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ABSTRACT 

To constrain the nature of the very first stars, we investigate the collapse and fragmentation of pri- 
mordial, metal-free gas clouds. We explore the physics of primordial star formation by means of three- 
dimensional simulations of the dark matter and gas components, using smoothed particle hydrodynamics, 
under a wide range of initial conditions, including the initial spin, the total mass of the halo, the redshift 
of virialization, the power spectrum of the DM fluctuations, the presence of HD cooling, and the number 
of particles employed in the simulation. We find characteristic values for the temperature, T a, few 
100 K, and the density, n ~ 10'^ — lO** cm~^, characterising the gas at the end of the initial free-fall 
phase. These values are rather insensitive to the initial conditions. The corresponding Jeans mass is 
Mj ~ 10'^ Mq. The existence of these characteristic values has a robust explanation in the microphysics 
of H2 cooling, connected to the minimum temperature that can be reached with the H2 coolant, and 
to the critical density at which the transition takes place between levels being populated according to 
NLTE, and according to LTE. 

In all cases, the gas dissipativcly settles into an irregular, central configuration which has a filamentary 
and knotty appearance. The fluid regions with the highest densities are the first to undergo runaway 
collapse due to gravitational instability, and to form clumps with initial masses ^ IQ^Mq, close to the 
characteristic Jeans scale. These results suggest that the first stars might have been quite massive, 
possibly even very massive with Mi, > 100 Mq. 

After a gas element has undergone runaway collapse, and has reached densities in excess of 10^ cm~'^, 
a sink particle is created. This procedure allows us to follow the evolution of the overall system beyond 
the point where the first nonlinear region would otherwise force the calculation to a halt. These later 
evolutionary stages, during which the clumps grow in mass due to accretion and merging with other 
clumps, are quite sensitive to the initial conditions. The key process in building up very massive clumps, 
with masses up to a few times W^Mq, is merging between clumps. Since the merging rate sensitively 
depends on the density of the gas, halos with the highest degree of central concentration are able to 
assemble the most massive clumps. Among these are halos with a low spin (A ~ 0.01), and with DM 
fluctuations imprinted according to a white-noise spectrum. 

Subject headings: cosmology: theory — early universe — galaxies: formation — stars: formation — 
hydrodynamics 



L INTRODUCTION 

This paper is concerned with the initial stages in the for- 
mation of cosmic structure. How did the universe evolve 
from the extreme uniformity of the big bang fireball into 
its highly organized and clustered present state? An in- 
creasing wealth of observational data has become available 
to guide theoretical thinking. The study of anisotropies 
in the cosmic microwave background (CMB) provides a 
window into the earliest phases of structure formation, 
when the primordial density fluctuations were still linear, 
and therefore amenable to an exact physical description 
(e.g., White, Scott, & Silk 1994; Lawrence, Scott, & White 
1999). Thus, we have a powerful probe of the last scatter- 
ing surface at z ~ 1000, corresponding to ~ 10^ yr after 
the big bang. Complementary to the CMB studies are ob- 
servations of high-redshift quasars and galaxies (e.g., Hu, 
Cowie, & McMahon 1998; Chen, Lanzetta, & Pascarelle 
1999; Fan et al. 2000). The light from these objects orig- 
inated at z ~ 5, when the universe was ~ 10^ yr old. 



Very little is known about the crucial era in between, at 
z — 1000 — 5, which has been termed the "dark ages" 
(e.g., Loeb 1999, Rees 1999). This serene epoch ended 
when the first luminous objects lit up the universe again. 
In the context of hierarchical scenarios of structure for- 
mation, as specified by a variant of the cold dark matter 
(CDM) model, the collapse of the first baryonic objects is 
expected at redshifts z ~ 50 — 10, involving dark matter 
halos of mass ~ IO^Mq (Tegmark et al. 1997). 

The importance of the first stars and quasars derives 
from the crucial feedback they exert on the intergalactic 
medium (IGM). A generation of stars which formed out of 
primordial, pure H/He gas (the so-called Population III) 
must have existed, since heavy elements can only be syn- 
thesized in the interior of stars. Population HI stars, then, 
were responsible for the initial enrichment of the IGM with 
heavy elements. From the absence of Gunn-Peterson ab- 
sorption in the spectra of high-redshift quasars, we know 
that the universe has undergone a reionization event at 
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z > 5 (Gunn & Peterson 1965). UV photons from the 
first stars, perhaps together with an early population of 
quasars, are expected to have contributed to the reion- 
ization of the IGM (Haiman & Loeb 1997; Ferrara 1998; 
Miralda-Escude, Haehnelt, & Rees 2000). The energy in- 
put from the first stars might have left a measurable im- 
print on the CMB on very small scales (Vishniac 1987; 
Dodelson & Jubas 1995). 

Despite an intense observational effort, the discovery 
of a true Population III star remains elusive. Surveys of 
the metal-poor population in the halo of our Galaxy have 
resulted in stars with metallicities Z ^ 10~^Zq (Beers 
2000). Spectroscopic studies of high-redshift Lyman a 
clouds, the supposedly most pristine objects in the uni- 
verse, find a minimum metallicity Z^^^ ^ 10~^.Zq (Cowie 
& Songaila 1998). Therefore, wherever we look, we find 
contaminated material. To probe the time when star for- 
mation first started, consequently entails observing at even 
higher redshift. This is the main purpose of the Next 
Generation Space Telescope (NGST) which is designed to 
reach ^ nJy sensitivity at near-infrared wavelengths (Loeb 
1998). In preparation for this upcoming observational rev- 
olution, the study of the first stars is very timely, provid- 
ing a theoretical framework for the interpretation of what 
NGST might discover, less than a decade from now. 

The question arises how one can make any progress in 
understanding primordial star formation, given the lack 
of direct observational constraints. The physics of the 
first stars, however, is characterized by some important 
simplifications, as compared to the extreme complexity 
of present-day star formation (Larson 1998; Loeb 1998). 
The absence of metals, and consequently of dust, leaves 
atomic and molecular hydrogen as the main agent of ra- 
diative cooling and the source of opacity. Magnetic fields 
were likely to be dynamically insignificant, prior to the 
onset of efficient (stellar) dynamo amplification (Kulsrud 
1997). The chemistry and heating of the primordial gas 
was not yet complicated by the presence of a UV radia- 
tion background. The IGM must have been a rather qui- 
escent place, with no source to sustain turbulent motion, 
as long as the density perturbations remain in their linear 
stage. Only after the explosion of the first supernovae, and 
the associated input of mechanical and thermal energy, is 
this state of primordial tranquility bound to change (Loeb 
& Haiman 1997; Miralda-Escude & Rees 1997). There- 
fore, the physics of primordial star formation is mainly 
governed by gravity, thermal pressure, and angular mo- 
mentum. This situation renders the problem theoretically 
more straightforward and tractable than the highly com- 
plex present-day case which continues to defy attempts 
to formulate a fundamental theory of star formation. Fi- 
nally, the initial conditions for the collapse of a primordial 
star forming cloud are given by the adopted model of cos- 
mological structure formation. The initial abundances of 
the chemical species are predicted by standard Big-Bang 
nucleosynthesis (e.g., Copi, Schramm, & Turner 1995). 

In this paper, we investigate the question: How do the 
first primordial star forming clouds evolve in the context 
of a hierarchical m,odel of structure form,ation? The col- 
lapse of the clouds, having masses close to the cosmological 
Jeans mass (~ lO^M©), results in the formation of high 
density clumps. These clumps are the immediate progen- 



itor of primordial protostars. This second stage in the 
primordial star formation process will be discussed in a 
subsequent paper (henceforth Paper II). 

The subject of the formation of the first stars has a 
long and venerable history (e.g., Yoneyama 1972; Hutchins 
1976; Silk 1977, 1983; Carlberg 1981; Kashlinsky & Rees 
1983; Palla, Salpeter, & Stabler 1983; Carr, Bond, & Ar- 
nett 1984; Couchman & Rees 1986; Haiman, Thoul, & 
Loeb 1996; Uehara et al. 1996; Tegmark et al. 1997; 
Omukai & Nishi 1998; Nakamura & Umemura 1999). Re- 
cently, it has become possible to address this problem in 
the context of full cosmological simulations, due to dra- 
matic improvements in numerical resolution, and in the 
modelling of the relevant gas physics (Anninos & Norman 
1996; Ostriker & Gnedin 1996; Abel et al. 1998; Abel, 
Bryan, & Norman 2000, ABN2000 henceforth; Fuller & 
Couchman 2000). 

Our approach is complementary to these last studies 
in that we focus on the collapse of an isolated overden- 
sity. This choice has its advantages and shortcomings. 
Paramount among the latter, we ignore the tidal effects 
of the large scale matter distribution which are responsi- 
ble for generating the angular momentum in the collapsing 
structures. The amount and distribution of angular mo- 
mentum therefore are free parameters in our simulations. 
In prescribing them, however, we can draw on the insight 
from cosmological numerical simulations. These result in 
a statistical description of the angular momentum (spin) 
which a given dark matter halo is expected to acquire (e.g., 
Barnes & Efstathiou 1987). The primary advantage of our 
method, on the other hand, is that it allows us to compre- 
hensively explore the behavior of the primordial gas under 
a variety of initial conditions. In doing so, we hope to 
single out the essential physics, and to test its robustness. 
In a previous publication, we have already presented first 
results (Bromm, Coppi, & Larson 1999). 

The organization of this paper is as follows. In §2, we 
discuss the relevant physical ingredients, both for the dark 
matter and baryonic components. A description of our 
numerical method is given in §3, whereas §4 presents the 
results of our exploratory survey. Finally, §5 contains our 
conclusions and debates avenues for further progress. 

2. THE PHYSICAL INGREDIENTS 

2.1. Dark Matter 

Current scenarios of cosmological structure formation 
assume the dark matter (DM) to be 'cold', in the sense 
of having negligible velocity dispersion. Candidates in- 
clude the lightest supersymmetric particle, the photino, of 
estimated mass ~ 100 GeV, and the invisil)le axion (e.g.. 
Peacock 1999). Primordial density fiuctuations in the cold 
dark matter have consequently survived on all scales due to 
the absence of free-streaming. Provided these fluctuations 
obey Gaussian statistics, as predicted by inflation, they 
are fully described by the power spectrum P{k), with k de- 
noting the wavenumber. The standard CDM model (e.g., 
Blumenthal et al. 1984) predicts that the fluctuations de- 
crease with mass, leading to a hierarchical (bottom-up) 
picture of structure formation. Variants of CDM agree 
qualitatively with each other, and we take the standard 
CDM model as a convenient template for our discussion. 
In the critical Einstein-de Sitter model, the rms density 
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contrast on a mass scale M grows in time according to 
a{M) = aQ{M) I [1 + z)^ as long as the fluctuation remains 
in the linear stage. fTo(M) is the linear extrapolation to 
the present epoch. In a model with a cosmological con- 
stant, growth is suppressed at recent epochs (z < 10), but 
at earlier times, the simple Einstein-dc Sitter behavior still 
applies. For a DM halo of given mass M, corresponding 
to a !/(T-peak in the Gaussian random field, the redshift of 
collapse can be estimated as 



1 + Zy 



(1) 



The threshold ovcrdcnsity for collapse is often taken to be 
5c = 1.69 (Peacock 1999). We discuss in Section 2.3 that 
the first stars are expected to form in DM halos of mass 
~ IQ^Mq, corresponding to 3 — 4cr peaks. At these small 
scales, halos of increasing mass collapse in rapid succes- 
sion. This ('cross-talk') behavior is characteristic of the 
CDM model (Rees 2000). The CDM power spectrum ap- 
proaches -P(fc) oc asymptotically on small scales. This 
is characteristic for the DM fiuctuations within a Popula- 
tion III star forming region. Since the baryonic mass is 
smaller than the initial Jeans mass, all perturbations in 
the gas have been wiped out by pressure forces. The pres- 
ence of the small-scale DM fluctuations, therefore, might 
play an important role in shaping the fate of the collapsing 
gas. 

A collapsing DM halo is expected to acquire angular 
momcntiim via tidal interactions with neighboring over- 
densities. The outcome of this process can be conveniently 
described by the dimensionless spin parameter 



A 



(2) 



where L, E, and M arc the total angular momentum, en- 
ergy, and mass, respectively, with G being Newton's con- 
stant. Numerical simulations result in an average spin 
parameter of A ~ 0.05 (e.g., Barnes & Efstathiou 1987). 

Once a perturbation on a given mass scale reaches over- 
densities 6 = {phaio — Pb) / Pb of close to Unity, with pf, being 
the density of the unperturbed background universe, the 
linear description breaks down, and the halo is undergo- 
ing nonlinear collapse. A convenient analytical model for 
the nonlinear evolution of DM halos is given by the spher- 
ical top- hat model (e.g., Padmanabhan 1993). At red- 
shifts close to Zyir, the DM particles reach a state of virial 
equilibrium. The approach to virial equilibrium occurs 
through the process often called violent relaxation in the 
literature (Lynden-Bell 1967). This process operates on a 
dynamical timescale, as opposed to the slow two-body re- 
laxation. The halo density after virialization is estimated 
to be 



(3) 



where po = 1-88 x 10~^g cm""* Vlmh is the density of 
the present-day universe. The baryons partake in the DM 
collapse, and acquire, through shocks and adiabatic com- 
pression, the virial temperature 
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Here, Rvir denotes the virial radius, fce is Boltzmann's 
constant, and mpj is the mass of a hydrogen atom. Al- 
though realistic DM halos are much more complicated 
than this simple model, the top-hat results give an intu- 
itive way to understand the physics of complex situations 
to within factors of a few. It also allows us to specify the 
initial conditions for our numerical simulations (see Sec- 
tion 4). 

In Figure 1, we show the density evolution of a top- hat 
perturbation, and compare the analytical prediction with 
the result of one of our numerical simulations. It can be 
seen that the prediction of pyir in equ. (3) is nicely repro- 
duced in the simulation, after an epoch of settling down 
to the equilibrium state. 
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Fig. 1. — Top-hat model for an overdcnsity collapsing at 
Zvir — 30. Shown is the number density of hydrogen atoms vs. red- 
shift. Solid line: Density evolution according to top-hat solution. 
Long-dashed line: Density evolution in the expanding background 
(E-dS) universe, pyir denotes the "virial plateau" , the estimated 
density after violent relaxation has established virial equilibrium. 
Diamond- shaped symbols: Evolution of the average gas density in 
adiabatic test case. The simulation is initialized at z = 100 to re- 
produce the top-hat collapse. As can be seen, the agreement with 
the analytical solution is good. Notice also that the average gas 
density settles eventually to a value close to pvir- 

2.2. Baryons 

The fate of the baryons is determined by the strife be- 
tween gravity and opposing pressure, and it is, therefore, 
essential to understand the thermal history of the g as, a 
topic which we address in the following subsection. 

2.2.1. Cooling and Heating 
The thermal evolution of the gas is governed by the 

F- A 

■-^^ + 



equation: 



~Dt 



P Dp 
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(5) 



(4) 



Here, P and p are the gas pressure and density, u is the spe- 
cific internal energy (in erg g""'^), and F and A are the con- 
tributions from radiative heating and cooling, respectively. 
The first term on the right-hand side describes adiabatic 
cooling due to expansion or heating due to compression. 
We now discuss the relevant processes of radiative cooling. 

In the absence of metals, H2 is the main coolant below 
~ 10'' K, which is the temperature range typically encoun- 
tered in collapsing Population III objects. Being homonu- 
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clear, H2 possesses no permanent dipole moment. Rota- 
tional transitions, therefore, can only occur via electric 
quadrupole radiation with the corresponding very small 
transition probabilities {A ~ 10~^^ s~^ for the lowest-lying 
transition). The H2 cooling fimction has been investigated 
by a number of authors whose results, however, differed by 
up to a factor of 100 (HoUenbach & McKee 1979; Lepp & 
Shull 1983). A crucial uncertainty concerned the low tem- 
perature (a few 100 K) regime, where quantum-mechanical 
effects were not yet taken into proper account. Recently, 
the cooling rates have been thoroughly recomputed. We 
have implemented H2 cooling using the new parameteriza- 
tion of Galli & Palla (1998). These authors have included 
improved collisional rates for T > 600 K (Martin, Schwarz, 
& Mandy 1996), as weU as for T < 600 K (Forrey et al. 
1997), and have assumed ortho- and para-H2 to be present 
in their equilibrium ratio of 3:1. In general, the cooling 
rate (in erg s~^ cm~^) can be written as 



A = ^ riuAuiAEui 



(6) 



u — 



where the sum extends over all possible transitions u I. 
Aui is the Einstein coefficient for spontaneous emission, 
AEui the respective energy difference, and n„ the density 
of H2 in the (upper) level u. In the low-density regime 
(n — > 0), each collisional excitation is almost instanta- 
neously followed by spontaneous emission. The level pop- 
ulations are then determined by the rate of H-H2 collisions: 



(7) 



with 7i„ being the coefficient of collisional excitation. At 
high density, the levels are populated according to LTE: 

n™(xnHexp(-A£;„,/fcBr)cxnH • (8) 
Summarizing the limiting behavior, one has: 



y^(LTE) ^ ^ ^ Ucrit 



(9) 



Bohr radius, oq — 0.5 A, one then finds for the lowest-lying 
(2 0) rotational transition the equivalent minimum tem- 
perature 



T = 



AE. 



20 



2mHao^s 



(12) 



resulting in Tmin ~ 500 K. The high-energy tail of the 
Maxwellian velocity distribution will allow the gas to at- 
tain somewhat lower temperatures, but it is one of the 
crucial aspects of the primordial gas physics that cooling 
down to temperatures below T ^ 100 K is not possible 
in the absence of coolants other than molecular hydrogen. 
The corresponding critical density for T ~ 100 — 300 K is 
then Ucrit — 10^ — 10* cm~^. These characteristic values 
of temperature and density, based on the microphysics of 
H2 cooling, leave their imprint on the thermal evolution of 
the primordial g will be discussed below. 




Fig. 2. — H2 cooling function (Galli & Palla 1998). Solid lines: 
Cooling rates per H2 molecule for various densities. From bottom 
to top: n = 10~^ , 10^ , 10^, 10^ cm""^. Diamond-shaped symbols: 
Cooling rate in LTE. This latter rate gives the maxiniuni possible 
cooling per molecule at a given temperature. Notice the saturation 
in the cooling rate due to the transition from NLTE to LTE level 
populations. 



The critical density ricrit, above which de-excitation due 
to collisions dominates over the competing radiative decay, 
is a function of temperature only, and is defined as 



f^crit 



A(LTE) 



(10) 
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The cooling function can finally be expressed in the form 



A 



A(LTE) 



1 + ricnt/nn 



(11) 



smoothly joining the low-density with the high-density 
case. In Figure 2, we show this cooling function for var- 
ious densities. At a given temperature, the cooling per 
molecule (A/riHa) H^st increases with density, but then 
approaches (for n > ricrit) the LTE value which consti- 
tutes the maximum possible cooling rate. Since H2 has a 
comparatively small moment of inertia, / ~ mnOQ, the ro- 
tational energy, Erot = L"^/^!, is correspondingly large for 
a given amount of angular momentum L = h^/J(J+l). 
Estimating the internuclear seperation to be of order the 



Another possible coolant in the primordial gas is HD, 
deuterium hydride. The low abundance, nn — 10~^nH, 
is partially offset by the fact that HD docs possess a per- 
manent electric dipole moment with the correspondingly 
larger radiative transition probabilities (A ^ 10~* s~^ for 
the lowest-lying rotational level, a factor ~ 1000 enhance- 
ment over H2). In addition, since the rotational transition 
1 ^ is allowed, HD cooling can reach lower tempera- 
tures than H2: Tmin — AEio/ks — 160 K. Consequently, 
cooling due to HD might become important at tempera- 
tures ^ 100 — 200 K. Our treatment of HD cooling relies 
on recent work of Flower et al. (2000), who have care- 
fully computed the relevant collisional rates, and provide 
an analytical fit to their results. 

Since radiative cooling to temperatures below that of 
the CMB 

TcMB = 2.7K(1 + z) (13) 

is thermodynamically not possible, we write for the cooling 
due to H2 and HD: 



A = A(T) - A(Tcmb) 



(14) 
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This approximate treatment ensures that T > Tcmb, un- 
less coohng proceeds via adiabatic expansion. 

Finally, we have included two cooling mechanisms of 
lesser importance. Firstly, for temperatures approaching 
^ 10^ K, cooling due to collisionally excited lines of atomic 
hydrogen becomes effective (Cen 1992): 

Ah = 7.5 X 10-^^(1+ T5/^)-iexp(-118,348/T)nen[H] 

(15) 

Only a small fraction of the gas in a Population III 
DM halo, however, reaches that high a temperature, and 
atomic hydrogen could not have facilitated the collapse 
and fragmentation of the primordial gas. Secondly, as 
long as there remains a residual degree of ionization, en- 
ergy is exchanged between the photons of the CMB and 
the electrons of the gas by (normal and inverse) Compton 
scattering: 

Acomp = 5.4 X 10-36(1 + ^) V(T - Tcmb) (16) 

At redshifts below z ^ 200, this coupling becomes increas- 
ingly weak, and does not significantly contribute to the 
cooling of the gas. 

In Figure 3, we show the relative importance of the 
various heating and cooling processes for a fluid element, 
whose time evolution is representative of what we find in 
our simulations. At early times, while the cloud is still 
expanding, cooling due to adiabatic expansion dominates. 
After turnaround, the gas heats up by adiabatic compres- 
sion. Only close to the moment of virialization docs H2 
cooling exceed the adiabatic heating, allowing the gas to 
cool upon further contraction. In these early stages, cool- 
ing due to HD is never important. After virialization, 
when the gas has settled into a cold (a few 100 K) cen- 
tral configuration, HD cooling becomes somewhat more 
important, without ever changing the thermal evolution 
dramatically. 
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Fig. 3. — Importance of cooling and heating terms during col- 
lapse. Solid lines: Cooling (in erg s^^ g^^) vs. time. Dotted 
lines: Heating (in erg s~^ g"^) vs. time. Time is measured in 
units of the initial free-fall time, ~ 2 X 10^ yr. Adiabatic cooling 
changes into adiabatic heating after turnaround. Compton heating 
(for T < Tcmb)i the other hand, turns into Compton cooling 
(for T > TQ]y[g). H2 is the dominant coolant, but can compete with 
the heating due to adiabatic compression only close to the moment 
of virialization. 

The treatment of radiative cooling, therefore, necessi- 
tates knowing the abundances of H2 and HD. This is the 



subject we take up next. 

2.2.2. Chemistry 

As shown in the previous section, radiative cooling is 
mainly due to H2. Cooling due to HD is less important, 
but we have included it for the sake of completeness. It 
is consequently essential to predict the respective abun- 
dances of these molecules. The appropriate primordial 
chemistry has the following simplifying features. Helium is 
assumed to be always completely neutral, and we neglect 
all reactions involving He, Hc+, Hc++. This assumption 
is justified since the temperature is typically low enough 
(T < 10* K) to render the He chemistry inert. In addi- 
tion, we also ignore all processes of photoionization and 
-destruction. Although photoreactions due to the CMB 
are very important at redshifts z > 100, they can be safely 
neglected in our simulations. Finally, there does not yet 
exist a UV background, prior to the onset of Population III 
star formation. Since H2 and HD arc formed in nonequi- 
librium, it is necessary to consider the respective reaction 
networks. First, we discuss the chemistry of H2 formation 
and destruction. 

Our H2 network is based on the compilation of Haiman 
et al. (1996), and is presented in Table 1. In the absence 
of dust grains, the main route for the production of H2 is 
given by the H" channel (McDowell 1961): 

R + e~ +hiy 

H" + H ^ H2 + e" 

An alternative formation channel relies on the intermedi- 
ary HJ. Since the H~ channel always dominates in our 
simulations, we ignore reactions involving H2'. The va- 
lidity of this assumption has been tested by comparison 
calculations with the full network. At redshifts z > 100, 
on the other hand, the HJ channel becomes important due 
to the ready destruction of H~ , which has a binding energy 
of only 0.75 eV, by energetic CMB photons. Thus, hydro- 
gen molecules are produced as long as there is a sufficient 
abundance of free electrons which act as catalysts in the 
H^ channel. For the low temperatures and weak shocks 
in the typical collapse of a primordial cloud, there is no 
efficient way to ionize the gas. The residual abundance of 
free electrons x = Ue/n ~ 10~* recombines on a timescale 
tree — {kreeXn)"^ ~ 10^ — 10*" yr. Thc recombination 
coefficient kree ^ 10^^^ cm'^ s^^ is evaluated for typical 
temperatures, T ~ 5000 K, and densities, n ~ 10^ cm"^, 
as found in our simulations. One can then estimate the 
maximal abundance of hydrogen molecules, / — nn^/n, 
according to 



The rate for the radiative attachment of H~ and e~ (see 
reaction (3) in Table 1) is approximately fcs ~ 10"-^^ cm^ 
s^^. Thc asymptotic abundance of H2, produced through 
the H^ channel, is then / ~ k^xntree — 10^"'. This predic- 
tion is borne out in our simulations under a wide range of 
circumstances. H2 abundances of 10"'' — 10"^, albeit low, 
nevertheless allow the primordial gas to efficiently cool. 
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thus enabling the condensation into high density struc- 
tures. Once the gas reaches densities n > 10*^ cm^'^, three- 
body reactions (Palla et al. 1983) convert the gas rapidly 
into fully molecular form. These reactions are included 
in Paper II, where we investigate the further collapse of a 
high density clump, but are not important here. 

The chemistry responsible for the formation and de- 
struction of HD has recently been discussed by Galli & 
Palla (1998). We adopt their minimum model, as listed in 
Table 2. 

2.3. Properties of the Star Forming Region 

To determine the properties of the primordial star form- 
ing region, two ingredients have to be considered. Firstly, 
one needs to know how the dark matter evolves. This his- 
tory of hierarchical gravitational clustering, resulting in 
the collapse of increasingly massive DM halos, is speci- 
fied by the adopted variant of the CDM model. Secondly, 
one has to address the question of whether the baryons 
are able to fall together with the dark matter. The op- 
posing effect of gas pressure in collapsing DM halos has 
been studied by Tegmark et al. (1997). These authors 
argue that the primordial gas can only undergo continued 
collapse and fragmentation, and consequently star forma- 
tion, if it manages to efficiently radiate away the gravi- 
tational binding energy which is released in this process. 
The efficiency of cooling can be quantified by the cooling 
timescale, tcooi — nk-QT / K, and is to be compared with 
the freefall time, f// — 1/y/Gp. One then has the classi- 
cal Rees-Ostriker criterion for fragmentation and contin- 
ued collapse (Recs & Ostrikcr 1977); tcooi < tff- In the 
early stages of the collapse, temperatures are low enough 
(~ 1000 K) for cooling to proceed mainly via the lowest- 
lying rotational transition of II2. In evaluating the crite- 
rion tcooi < iff, Tegmark ct al. (1997) find the minimum 
virial temperature necessary for sufficient cooling, and by 
applying equ. (4) the corresponding minimum halo mass 
as a function of collapse redshift. The result of this calcu- 
lation is that a 3(j-peak of mass ~ 10^ Mq, and collapsing 
at a redshift Zyir ~ 30, does meet the requirement for effi- 
cient cooling. Halos with these characteristic masses and 
collapse redshifts are therefore predicted to be the sites for 
the formation of the first stars. 

To understand the complicated physics of gas fragmen- 
tation, we now turn to a description of our numerical sim- 
ulations. 

3. NUMERICAL METHODOLOGY 

In this section, wc describe the elements of our nu- 
merical approach. The evolution of the dark matter and 
gas components is calculated with a version of TREESPH 
(Hernquist & Katz 1989), combining the smoothed particle 
hydrodynamics (SPH) method with a hierarchical (tree) 
gravity solver. The details of the gravitational N-body 
solver and the treatment of the hydrodynamics are dis- 
cussed in the Appendix. We here present our additions 
to the code which are necessary for the investigation of 
the primordial gas. These are a method to solve the pri- 
mordial chemistry network, and a technique to create sink 
particles. 

3.1. Solving the Reaction Network 



Following the abundance evolution of the 8 species H, 
11+ , H", H2, D, D+, HD, and c", entails solving the 
corresponding coupled set of kinetic equations. Due to 
the very short reaction timescales, compared to the dy- 
namical time, the adopted method of solution has to be 
implicit, as in the case of the thermal energy equation. 
Traditional matrix-based techniques such as the STIFBS 
routine (Press et al. 1992), or the Livermore LSODAR 
solver (Hindmarsh 1983), prove to be computationally too 
expensive for three-dimensional applications. We there- 
fore have implemented a fast method of solution which is 
based on the approximate backwards differencing formula 
(BDF). This approach has been pioneered by Anninos et 
al. (1997). Consider the rate equation for a given species 
i, which can be expressed as 



dt 



—1 = -Dm + C 



(17) 



All reactions contributing to the destruction of species i 
are contained in D, while all reactions leading to the pro- 
duction of that species are summarized by the symbol C. 
Both D and C are functions of temperature and the abun- 
dances of the other species. Evaluating the right-hand 
side of equation (17) at the new timestep, results in the 
discretization 



^old 



1 -I- £)"™Ai 



(18) 



We update the hydrogen and deuterium species in the fol- 
lowing order, which has been found by experimentation 
to be both stable and accurate: e^, H+, H^, H2, H, D+, 
HD, D. An estimate for C"'="' and is obtained by 

first inserting the abundances at the old timestep, and 
then successively replacing them with the updated ones, 
as the algorithm proceeds from species to species. In de- 
termining the abundances of the normal hydrogen species, 
we do not include the reactions of the deuterium network. 
This is justified by the very low abundance of deuterium 
(riD ^ 10~^nH)- The evolution of the deuterium network, 
on the other hand, crucially depends on the abundances 
of the normal hydrogen species. 

For increased accuracy, a subcycling procedure is 
adopted: The algorithm loops over the species repeatedly, 
each time with a small timestep Ats„6 given by 



(19) 



where n,. is the abundance of free electrons, and e = 0.1 is 
an accuracy parameter. 

The reactions determining the abundance of the inter- 
mediary species H~ proceed much more rapidly than the 
other reactions. Therefore, one can safely assume that H~ 
is always present with its equilibrium value: 



(^4 + A;ii)nH + (fcs + A;i2)njj+ -|- kioUe 



(20) 



The identification of the reaction rates can be found in 
Table 1. 

To test our chemistry solver, we compare the analyti- 
cally determined equilibrium abundances as a function of 
temperature, obtained by setting dui/dt = for all i, with 
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those given by the chemistry solver, after the integration 
has converged. The analytical estimates are derived as 
follows. For H"'", we assume Ue = n 
coUisional balance: 



JJ+ and consider the 



H + e" 
H+ + e- 



H+ + 2e- 



With the reaction rates in Table 1, one finds 



ki 



njj+ = — TiH 



(21) 



The abundance of atomic hydrogen is assumed to be 
nn = nu,tot — '^h+' where nu,tot is the total density of 
normal hydrogen. For H2, we obtain (see Table 1): 



nii2 



(fce + fcs + k9)n-a + fcyriH 



(22) 



The abundance of the intermediary species is given by 
equation (20), evaluated with riua — 0. For the deuterium 
species, D+ and HD, the corresponding expressions are: 



and 



fcj)2^D^H+ 

koirie + kDsnu 



(23) 



(24) 



Here, koi^ko^, ■ ■ ■ correspond to the reactions in Table 
2, and we assume no = nu,tot — ''^0+' with nu,tot being 
the total density of heavy hydrogen. The result of this 
comparison is presented in Figure 4. As can be seen, the 
chemistry solver nicely reproduces the analytical expecta- 
tion. 



1.0 
0.8 
0.6 
0.4 
0.2 
0.0 



1.0 
0.8 
0.6 
0.4 
0.2 
0.0 



T [K] 



T [K] 

HD 













'1 


< 





1000 



T [K] 



1000 



T [K] 



Fig. 4. — Testing the chemistry solver. Shown are the equiUb- 
rium abundances for the species H+, H2, D+, and HD. In each case, 
fractional abundance is plotted vs. temperature. Solid lines: Abun- 
dances from allowing the chemistry solver to reach convergence (to 
within 10"^"). Diamond- shaped symbols: Analytic estimate for the 
equilibrium value. It can be seen that the equilibrium abundances 
are nicely reproduced by the solver. Notice also the extreme simi- 
larity between the normal hydrogen and the deuterium species. 



3.2. Creation of Sink Particles 

When in the course of a simulation the gas attains in- 
creasingly high density, the SPH smoothing length de- 
creases according to /i ~ ^neigh''^~^^^- "^^^ Courant con- 
dition then enforces the adoption of smaller and smaller 
tinicsteps At ~ h/cs, with being the somid speed. Upon 
the onset of gravitational instability and the resulting run- 
away collapse, the simulation therefore effectively grinds 
to a halt. To overcome this fundamental numerical limita- 
tion, and to be able to follow the evolution of the overall 
system for a few dynamical times, a method of creating 
sink particles is necessary. Here, we describe our strategy 
in doing so. Recently, Bate, Bonnell, and Price (1995) 
have developed another such technique in the context of 
SPH and applied it to the study of protobinary accretion. 

In developing the algorithm, one has to address the fol- 
lowing questions. (1) When a sink particle is created dur- 
ing a simulation, would the incorporated particles really 
continue to collapse, or would they escape again from each 
other, if one were to follow their evolution further? To en- 
sure that only gravitationally bound particles are merged 
to form a sink particle, we utilize the runaway nature of 
the Jeans instability. Consider how the density does evolve 
with time for the SPH particles close to the first collapsing 
region. Densities are initially close to n ^ 10* — 10^ cm""^, 
and subsequently either experience only a modest rise in 
density, or a rapid increase by many orders of magnitude. 
This dichotomy allows for the unambiguous identification 
of the collapsing gas particles. The first, and physically 
most important, criterion for a particle to be eligible for 
merging therefore is: (a) n > nth- The value of nth has to 
be adjusted to the physical characteristics of the problem. 
We find that nth — 10* cm~^ works well in our case of pri- 
mordial gas which collapses and fragments in the center 
of a dark matter halo. Another reason for this choice lies 
in the fact that beyond a density of 10^ cm~^, additional 
physics has to be considered: Three-body reactions con- 
vert the gas into fully molecular form, and opacity effects 
start to become important. We address this high-density 
regime in Paper II, but otherwise assume that n < 10* 
cm~^. A second test a particle has to pass is: (b) V-v < 0, 
i.e., that it is part of a converging flow. 

To further examine whether a given collection of parti- 
cles is gravitationally bound, one has to determine whether 
the binding energy of the system is negative (criterion (c)): 



Etotal — Eqrav + Ekin + Eth < 



(25) 



where Egrav, Ekin, and Eth are the total gravitational, 
kinetic, and thermal energies, respectively. A system 
will continue to collapse only if gravity overwhelms the 
combined opposing efforts of thermal pressure and rota- 
tion. Defining the usual parameters a = Eth/\Egrav\ and 
/3 = Erot/\Egrav\, With Erot being the total rotational en- 
ergy of the system, the requirement of continued collapse 
can be written as (criterion (d)): 



a + (3<l 



(26) 



In Figure 5, we show these quantities as a function of dis- 
tance from the density maximum, at the instant of creating 
the sink particle. 
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Fig. 5. — Search radius for inclusion into sink particle. Solid 
line: Total energy vs. distance from density maximum. The 
total energy of all particles within a given radius is defined as 
^total = Egrav + Ekin + -^t/t > a^id plotted with an arbitrary normal- 
ization. Distances are in units of h(rmax), the smoothing length of 
the particle with the highest density. Dotted line: Sum of the ratios 
a = Eth/\Egrav\ and /3 = Erot/\Egrav\ vs. r/h. Only particles 
within Ar < 2h {dashed line) are considered to become part of the 
sink particle. It is evident that this is a very conservative choice, 
safely meeting the criteria a+/3 < 1 and E^^^i^i < 0. A more aggres- 
sive merging would search within Ar < lOh, but this would entail a 
more elaborate testing of whether a candidate particle would really 
end up in the clump, or escape again. 



rci 



m,;r,: 



Mci 

and a mass-weighted velocity 



vci 



Mci 



(28) 



(29) 



where the summation extends over all contributing par- 
ticles. Wc treat the sink particle as being still gaseous, 
and as consequently still participating in the SPH inter- 
actions and smoothing procedure. A constant gas den- 
sity, nci = nth, and temperature, Tci = J^i'^i'^i/Mci, 
are assigned to the clump. The resulting (constant) pres- 
sure is higher than in the surroundings which is slightly 
less dense and hot. This prescription approximately mod- 
els a negative pressure gradient. The unphysical infall 
of neighboring particles onto the clump is therefore pre- 
vented. Such a situation would occur if the sink parti- 
cle were assumed to be coUisionless, unless special care is 
taken to formulate appropriate boundary conditions (see 
Bate et al. (1995)). With continuing mass accretion, grav- 
ity becomes increasingly dominant, and the sink particle 
asymptotically assumes a fully coUisionless character. In 
Figure 6, we demonstrate that this procedure leads to a 
satisfactory treatment of the boundary between the clump 
and its neighboring particles. 



Our merging algorithm now proceeds by sorting the eli- 
gible particles according to density. Subsequently, a search 
is performed within a range of r search = '^hmax around the 
particle with the highest density. All particles within this 
volume are merged to form the sink particle, provided they 
fulfill criteria (a) and (b). This procedure is repeated un- 
til all eligible particles are assigned to a sink particle. Wc 
assume that passing criteria (a) and (b) implies Efofai < 
and a + /3 < 1, without explicitly testing for it. This as- 
sumption is justified by the following physical argument. 
As can be seen in Figure 5, the choice of r search = '2hmax is 
very conservative. All particles within this search volume 
are safely bound, and far from pressure or rotational sup- 
port. Surrounding the search volume, there is a massive, 
collapsing envelope which is part of the Jeans unstable 
flow. Consequently, a more aggressive merging could be 
adopted with r search — Whmax- To reliably determine 
whether one of the less bound, outer particles is to be 
merged, the explicit testing for criteria (c) and (d) would 
be necessary in this case (Bate et al. 1995). We have 
carefully verified the reliability of our merging procedure 
by performing test calculations with and without criteria 
(c) and (d), and find that in each case our assumption, 
(a)-h(b)^ (c)-h(d), is valid. 

The next conceptual step concerns the properties of the 
newly created sink particle, and how to treat the bound- 
ary to the neighboring region. This leads to the second 
question: (2) Does the presence of a sink particle unphys- 
ically affect the SPH particles in its vicinity? Upon its 
formation, the sink particle (or clump) has mass 

Mci=^mi , (27) 
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Fig. 6. — Effect of merging on neighboring particles. Right 
column: Properties of SPH particles in the vicinity of the sink par- 
ticle, briefly after its creation. Left column: Particle properties for 
the comparison case, where no merging takes place, and the gas 
evolution is followed to increasingly higher density. Top panels: 
Number density vs. distance from density maximum. Distances arc 
in units of h(rmax), the smoothing length of the particle with the 
largest density, evaluated at the instant of merging. Middle pan- 
els: Smoothing length vs. r/h. Bottom panels: Gas pressure vs. 
r/h. The vertical lines at r/h = 2 delineate the effective size of the 
sink particle. It can be seen that the thermodynamic properties, 
together with the radial gradients of density and pressure, are very 
similar in the two cases. The smoothing lengths, on the other hand, 
are noticeably larger in the presence of a sink particle. This is to 
be expected, since removal of particles entails a degrading spatial 
resolution. 

The final element of the algorithm handles the subse- 
quent accretion of gas, and addresses the question: (3) 
Is the accretion of a given gas particle physically justi- 
fied, or does the particle venture only temporarily into the 
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vicinity of the clump! A particle is accreted onto a pre- 
existing clump if criteria (a) and (b) are fulfilled, and if 
it approaches to within the accretion radius ta = ^hci, 
where hci is the smoothing length of the accreting clump. 
This procedure proves to be very reliable, again due to 
the fact that in order to reach a density of 10^ cm"'^, the 
particle has to be part of a Jeans unstable gas flow which 
has already proceeded well into its runaway collapse. The 
position and velocity of the merged clump and accreted 
particle are again taken to be the mass- weighted averages. 
The density and temperature of the clump retain their 
initial values, as mentioned above. 

This algorithm allows us to investigate the complex dy- 
namics of clump formation, subsequent gas accretion, and 
the occasional merger of clumps. 

4. THE SIMULATIONS 

With all the ingredients in hand, we now turn to the 
description of our simulations, and present an exploratory 
survey of how the primordial gas behaves under a wide 
range of initial conditions. We begin by describing our 
procedure of initializing the numerical simulations, and 
explaining our choice of parameters. We then turn to the 
description of the simulations, focusing first on the early 
evolutionary stages, and subsequently on the later ones. 

4.1. Initial Conditions 

We model the site of primordial star formation as an 
isolated overdensity, corresponding to a high— cr peak in 
the random field of density perturbations. The numeri- 
cal simulations are initialized at Zi = 100 by endowing a 
spherical region, containing dark matter and gas, with a 
uniform density and a Hubble expansion according to the 
top-hat solution (see Section 2.1). The top- hat overden- 
sity is embedded in an Einstein-deSitter universe with a 
Hubble constant h = Fo/(100km/ s / Mpc) = 0.5. We 
consider halos collapsing at Zyir ~ 20 and at z^r ~ 30, 
corresponding to ~ 2a and ~ 3a peaks, respectively, with 
total masses of 2 x lO^M© and 2 x IO^Mq. The baryonic 
mass fraction is usually taken to be Ob = 0.05, but we 
also consider the case = 0.20. 

While this is clearly a rather idealized initial configura- 
tion, the subsequent evolution of the halo quickly departs 
from a simple, monolithic collapse. In response to the im- 
printed small-scale DM fluctuations (see below), the halo 
develops a very lumpy, nonspherical morphology. This 
model, which allows for the formation of DM substructure, 
and the hierarchical merger of DM clumps, embedded in 
a collapsing background medium, does qualitatively de- 
scribe the dynamical behavior of the cold dark matter on 
very small scales. 

The major shortcoming of this approach lies in the in- 
ability to self-consistently incorporate the angular momen- 
tum of the halo. In a realistic cosmological setting, angular 
momentum is generated by the tidal interaction of neigh- 
boring overdensities which are, in general, nonspherical. 
For an isolated halo, one therefore has to specify the distri- 
bution and magnitude of the angular momentum in an ex- 
plicit way. We here assume that the halo is initially in rigid 
rotation with a given angular velocity w. We prescribe lo 
in accordance with the prediction for the spin parameter 
A, as found in cosmological N-body simulations (see Sec- 
tion 2.1). In computational units with G = M = R = 1, 



one approximately has the relation: lu ~ 6.7A. The initial 
angular velocity has values uj = 0.1,0.2, and 0.4, corre- 
sponding to A = 0.015,0.03, and 0.06, respectively. 

We now describe our procedure of initializing the posi- 
tions and velocities of the DM and gas particles. To im- 
print small-scale density fluctuations on the dark matter, 
we carry out the following steps. After setting up the DM 
particles on a Cartesian grid, we impose periodic bound- 
ary conditions, and perturb the particles according to the 
Zcldovich (1970) approximation. The random density field 
S{x) has the Fourier decomposition S = Y^Sj: cxp{ik ■ x) 

with Sj: = ^^exp(i(/?g). A given mode with wavevector k is 
specified by the random phase (p^, distributed uniformly 
in the interval [0, 27r], and the amplitude = v^Pik), 
where v is drawn from a Rayleigh distribution (see Pad- 
manabhan 1993). The particles, having imdisturbed po- 
sition q, are displaced in such a way as to reproduce 
the desired density field. In the Zeldovich approach, one 
finds for the displaced particle position: x = q + f{q)- 
The displacement field / is related to the density via 
5 = — V • /. It is then straightforward to show that 
/ = X^/j:exp(ifc • q) with /g = i5j:k/k^. The Zeldovich 
approximation also allows one to self-consistently assign 
peculiar velocities: Vpe.c.i = Hif ^ where Hi is the Hub- 
ble parameter at z,; = 100. Adding the peculiar veloc- 
ity to the Hubble expansion and the solid-body rotation, 
gives the resulting initial velocity for each DM particle: 

Vi = VH,i + Vrot,i + Vpec,i- 

For Gaussian fluctuations, the power spectrum P{k) = 
|(5^P = Ak" fully describes the random density held. On 
small scales (M < lO^M©), the standard CDM model pre- 
dicts an asymptotic behavior of P{k) cx k^"^, and we take 
n = —3 as our standard value. To investigate the depen- 
dence of the gas fragmentation on the character of the DM 
substructure, we also consider the case P{k) oc fc", corre- 
sponding to white-noise perturbations. To flnally fix the 
amplitude A, we specify the initial variance of the fluctu- 
ations 

af=A^k^ . (30) 

The summation is over all contributing modes, where the 
minimum wavenumber is given by the overall size of the 
Cartesian box, and kmax by the Nyquist frequency. We 
typically have af ~ 0.1. 

The rms fluctuation at the moment of collapse is then 

.(. = 30).(l±il)...l . 

This choice ensures that the substructure develops on a 
similar timescale as the overall collapse of the background 
medium. A potential problem with this initialization pro- 
cedure is the poor sampling of k-space for the longest 
wavelength modes. Since the overall structure of the halo 
is predominantly determined by these modes (for a k~^ 
spectrum), the DM morphology will change significantly 
from one random realization to the next, and we can there- 
fore not claim to have simulated the DM morphology in a 
representative way. This is true in particular for the early 
stages of the collapse, whereas at later times, different ini- 
tial configurations approach similar equilibrium states. To 
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remedy this shortcoming, we have studied cases with dif- 
ferent random rcahzations of the DM fluctuations, keeping 
all other parameters constant. Finally, all particles within 
a given radius are selected for the simulation. The typical 
number of DM particles is A'dm — 17, 000, but we have 
also performed test calculations with N-qm — 130, 000. 

The SPH particles, representing the baryonic compo- 
nent, are placed randomly within the given spherical vol- 
ume. This procedure inevitably introduces unphysical 
shot-noise. As opposed to the dark matter, however, these 
numerically induced perturbations are quickly erased by 
pressure forces, since the initial gas mass is close to the 
Jeans mass, Mj ~ lO^M©. The gas particles have ve- 
locities, consisting of Hubble expansion and rigid rota- 
tion: Vi = VH.i + Vroi.i- The fractional free electron and 
hydrogen molecule abundances are taken to be, respec- 
tively, Xi = 4.6 X 10"^ and = 2 x 10"^ (Anninos & 
Norman 1996). We assume a deuterium abundance of 
no = 4 X 10~^nH (Galh & Palla 1998), and initiahze the 
density of D+ and HD according to 71^+ = 10~^nH, and 
?^HD = 10~^nH2j respectively. The gas temperature is fi- 
nally chosen to be Tj = 200 K. The typical number of SPH 
particles is A"sph — 16, 384, but we have also performed 
simulations with iVspH — 65,536 and Nsph — 131,072, 
to test for convergence. 

Initializing our simulations at Zi = 100 poses the ques- 
tion whether this is early enough to investigate objects 
that are collapsing at z ^ 30. The crucial initial values for 
the chemical abundances, the gas temperature and density 
are actually obtained by integrating the respective gov- 
erning equations beginning at 2; = 1100, along the line of 
Tegmark et al. (1997). Although the use of the Zeldovich 
approximation in advancing the density perturbations to 
Zi = 100 will introduce an error, the overall nature of the 
DM collapse would not change in a significant way if the 
simulations were started at somewhat higher redshift. 

In Table 3, we summarize the parameters of the differ- 
ent simiilations. Run A constitutes the fiducial case, which 
we will describe first, and against which we subsequently 
compare the other runs. 

4.2. Early Evolution 

Given the almost complete absence of observations to 
constrain the theoretical study of primordial star forma- 
tion, it is important to single out the essential physical 
processes. We here present evidence for the existence of 
characteristic values for the temperature and density of 
the primordial gas, which in turn translate into a charac- 
teristic Jeans scale for fragmentation. To illustrate this 
characteristic behavior, we first describe a fiducial simu- 
lation (Run A) in greater detail, and then investigate the 
effect of varying the initial conditions. 

4.2.1. The Characteristic Mass Scale for Fragmentation 

The parameters of Run A, a halo of total mass 2 x lO^M© 
with IO^Mq in baryons, and virializing at Zyir — 30, are 
chosen to safely satisfy the Rees-Ostriker criterion for con- 
tinued collapse and fragmentation, tcooi < tff (sec Section 
2.3). In a previous publication, we have described a simi- 
lar case (Bromm, Coppi, & Larson 1999). The simulation 
presented here has improved on that earlier treatment by 
including cooling due to HD, and a more realistic initial 



H2 abundance of = 2 x 10^® (instead of fi = 10^ ). 
Furthermore, Run A is initialized with a different random 
realization of the DM fluctuation field. 
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Fig. 7. Run A: Initial configuration for top-hat collapse. The 
halo has a total mass of 2 x lO^^Afrr, , and is endowed with a Hub- 
ble expansion such that virialization occurs at Zyir — 30. Top tow: 
The DM particles are perturbed from a regular grid according to 
P{k) oc k~^. Bottom row: The gas particles are placed at random, 
and comprise a mass fraction of Q^b = 0.05. Both components are 
initially in solid body rotation with the angular momentum vector 
pointing in the ^^-direction. Left panels: Face-on view. Right panels: 
Edge-on view. 
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Fig. 8. — Run A: Morphology at z = 33.5. The manner of 
presentation is the same as in Figure 7. The DM has developed 
significant substructure, and the baryons are just beginning to fall 
into the corresponding potential wells. 
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Fig. 9. — Run A: Morphology at z = 31.2. The convention in 
Fig. 7 is adopted for the rows and columns. The box size is 30 pc. 
The DM is in the process of undergoing violent relaxation with the 
concurrent smoothing out of the substructure. Having developed a 
lumpy and elongated morphology, the gas has settled into the center 
of the DM potential well. Shown is the situation briefly after the 
formation of the first clump of mass 1400AfQ {small box). 
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Fig. 10. — Run A: Gas properties at 2 = 31.2. (a) Free electron 
abundance vs. hydrogen number density (in cm~^). At densities 
exceeding n ^ 10^ cm~^, recombination is very efficient, and the 
gas becomes almost neutral, (b) Hydrogen molecule abundance vs. 
number density. After a quick initial rise, the H2 abundance ap- 
proaches the asymptotic value of / r^ 10~^, due to the operation 
of the H~ channel, (c) Gas temperature vs. number density. At 
densities below ~ 1 cm~^ , the gas temperature rises because of adi- 
abatic compression until it reaches the virial value of T„j^ ~ 5000 K. 
At higher densities, cooling due to H2 drives the temperature down 
again, until the gas settles into a quasi- hydrostatic state at T r^ 500 
K and n ~ 10'* cm~^. Upon further compression due to the onset of 
the gravitational instability, the temperature experiences a modest 
rise again, (d) Jeans mass (in Mq) vs. number density. The Jeans 
mass reaches a value of Mj r^J 10^ M© for the quasi-hydrostatic gas 
in the center of the DM potential well, and reaches the resolution 
limit of the simulation, Mres ~ 2OOM0 , for densities close to the 
merging threshold of n = 10*^ cm^''. 

In Figure 7, we show the initial configuration of Run 
A at Zi = 100. Comparing the dark matter and gaseous 



components, the different way of initializing them, as de- 
scribed above, is clearly visible. Initially, the halo is still 
expanding, until the moment of turnaround at Zta — 50. 
The subsequent collapse of the dark matter, ocurring on a 
dynamical timescale of tj-j ^ 1/ ^Gp(z = 100) 5 x lO''' 
yr, leads to the eventual establishment of virial equilib- 
rium at Zvir ^ 30. Since the initial gas pressure is dynam- 
ically unimportant, the baryons freely fall together with 
the dark matter. Upon compression, the gas temperature 
rises adiabatically until it reaches the virial value 



GMmn 

2/cb Rvir 



5000 K 



(31) 



where Ryir — 100 pc is the virial radius. At this point, 

enough H2 molecules have been formed (/ ^ a few 10~^) 
to provide an efhcient cooling mechanism. Consequently, 
the temperature decreases again with further compression. 
Figure 8 shows the situation at z = 33.5, briefly before 
the virialization of the dark matter. In response to the 
initially imprinted fc~'^-noisc, the dark matter has devel- 
oped a pronounced substructure. Evidently, the collapse 
proceeds in a very inhomogeneous manner, far from the 
monolithic, spherically symmetric evolution of the analytic 
top-hat model. The baryons have just begun to fall into 
the potential wells which are created by the DM substruc- 
ture. Thus, the DM imparts a 'gravitational head-start' 
to certain regions of the gas, which subsequently act as 
the nucleization centers for the formation of high-density 
clumps. 

At the end of the free-fall phase, shown in Figure 9, the 
gas has developed a very lumpy, filamentary structure in 
the center of the DM potential. By now, the dark matter 
has lost the memory of the primordial perturbations, but 
only after having imprinted its signature on the gas. This 
almost complete erasure of the DM substructure might be 
due to insufficient numerical resolution, as has been re- 
cently suggested by, e.g., Moore et al. (1999). We have 
performed a test calculation with eight times as many DM 
particles (A'dm ^ 130,000 instead of 17,000), and find no 
qualitative change, although a still larger number of parti- 
cles may be required (10^ — 10^) to see the survival of the 
DM lumps. 

The corresponding thermodynamic and chemical state 
of the gas is summarized in Figure 10. Since the abun- 
dances, temperature and density are plotted for every SPH 
particle, this mode of presentation has an additional di- 
mension of information: Particles accumulate ('pile up') 
in those regions of the diagram where the evolutionary 
timescale is slow. In panel (c) of Figure 10, one can clearly 
discern such a preferred state at temperatures of a few 
100 K, and densities of 10^ — 10^ cm^^. These character- 
istic values have a straightforward physical explanation, 
along the line of argument presented in Section 2.2 re- 
garding the microphysics of II2 cooling. A temperature of 
T ~ 100 — 200 K is the minimum one attainable via H2 
cooling. The corresponding critical density, beyond which 
the II2 rotational levels are populated according to LTE, 
is then ricrit — 10"^ — 10^ cm~^. Due to the now ineffi- 
cient cooling, the gas 'loiters' and passes through a phase 
of quasi-hydrostatic, slow contraction. 
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Fig. 11. — Run A: Important timcscalcs at z = 31.3. Shown 
is the situation briefly before the formation of the first clump, (a) 
Free-fall timescale {solid line) and cooling timescale {dotted sym- 
bols) vs. number density (in cm"^). Timescales are in units of 10^ 
yr. The gas particles pile up at a density of n ~ 10^ cm~^, where 
^cool —iff- (t") Free-fall timescale {solid line) and sound-crossing 
timescale {dotted symbols) vs. number density (in cm^"^). The on- 
set of the Jeans instability (i.e., tsound > *//) " ~ 10^ cm~^ 
coincides with the condition tcool — iff in panel (a) . 



Further insight can be gained by considering the char- 
acteristic timescales of the problem, which are displayed 
in Figure 11. We consider the free-fall time tff, the cool- 
ing time tcool — nkBT/A}i2, and the sound-crossing time 
tsound — Lchar/cs- Hcrc, Lchar ~ 1 pc is thc characteristic 
size of the filamentary gas in Figure 9. The comparison 
of tcool and tff explains the evolutionary behavior of the 
gas, as described above. At densities exceeding n ~ 10° 
cm~^, the gas can cool efficiently (i.e., tcooi < tff), until it 
reaches the quasi- hydrostatic phase, where tcooi ^ tff. To 
move away from this loitering regime, and to attain higher 
densities, the gas has to become gravitationally unstable. 
The condition for the onset of instability, t sound > tff, is 
shown in panel (b) of Figure 11. The gas is Jeans unsta- 
ble for n > 10* cm^'^. Alternatively, we can evaluate the 
Jeans mass for the characteristic values T ~ 200 K and 
n ~ 10^ - 10* cm-3, resulting in Mj ~ IO^Mq. When 
enough gas has accumulated in a given region to satisfy 
M > Mj, runaway collapse of that fluid region ensues. We 
find that the gas becomes self-gravitating (pe > Pdm) co- 
incident with the onset of the Jeans instability. Although 
the dark matter has played an important role in determin- 
ing where most of the gas ends up, it henceforth ceases to 
influence the primordial gas on its further course to star- 
dom. 



4.2.2. The Onset of Gravitational Instability 

After the onset of instability, the temperature rises 
again, up to T ~ 1000 K at n ~ 10^ cm~^, but not suf- 
ficiently to halt the collapse. Once the gas has reached a 
density of nth ~ 10^ cm~^, a clump (or sink particle) is 
formed, and we do not follow the evolution to increasingly 
higher density. The initial mass of the clump, as shown 
in Figure 9, is 14OOM0, close to the characteristic Jeans 
mass of Mj ~ lO^M©. 
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Fig. 12. — History of first runaway fluid element, (a) Gas tem- 
perature (in K) vs. cosmic time (in 10^ yr). Dashed line: Charac- 
teristic temperature Tchar — 300 K. (b) Hydrogen number density 
(in cm^"^) vs. cosmic time. Dashed line: Characteristic density 
f^char — 10^ cm^'^. Tlic fluid clement spends ~ 10^' yr at temper- 
atures and densities close to the characteristic values. After this 
period of "loitering", the runaway collapse sets in, operating on a 
timescale ~ 10^ yr. 

In Paper II, we take up the question of how the fur- 
ther collapse of a clump does proceed, up to densities of 
n ^ 10^^ cm^''. In doing so, we find no indication for 
further subfragmentation. Despite the additional boost 
in the cooling due to the action of three-body reactions, 
which convert the gas into almost fully molecular form at 
n > 10^ cm^"^, no runaway cooling occurs. At the end 
of the simulation presented in Paper II, a central core of 
~ IOOM0 is in a state of free-fall, surrounded by an ex- 
tended envelope with an approximately isothermal density 
profile, p oc r~^. 

Another way to understand the nature of the runaway 
collapse is shown in Figure 12, where we plot the history 
of the first runaway SPH particle. This particle marks the 
center of the fluid region which is first to become grav- 
itationally unstable. It is evident that the temperature 
reaches T ~ 300 K, and subsequently oscillates around 
that value for approximately 10^ yr, to finally experience 
a very rapid rise to ^ 1000 K. The oscillatory behavior 
is due to the negative feedback of the H2 cooling, with 
increased cooling for higher temperature, and vice versa. 
Similarly, the density stays at n ~ 10** cm~'^ for ~ 10^ yr, 
before the onset of the runaway collapse. 

It has been pointed out by Bate & Burkert (1997) that 
to avoid numerical fragmentation, one has to resolve the 
Jeans scale, Mj > Mres- One can estimate the resolution 
limit, Mres, of the simulation as 

Mres^f^^We . (32) 

For Run A, where A^sph = 16, 384, this results in Mj-es ^ 
2OOM0, and the criterion above is reasonably well satis- 
fied. To test for numerical convergence, we have performed 
a simulation with A^sph = 131,072 (Run B), and a cor- 
responding mass resolution of Mres ~ 2^Mq. This high- 
resolution run confirms that clumps initially form with 
masses close to IO^Mq. 
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It is a longstanding question whether the Jeans mass is 
relevant for the understanding of the characteristic mass of 
present-day star formation. In the primordial case, how- 
ever, where one can argue that magnetic and turbulent 
pressures are initially unimportant, one is left with the 
classic battle between gravity and thermal pressure, as 
originally envisioned by Jeans (1902). Primordial star for- 
mation, therefore, might be the most clear-cut setting for 
the application of the Jeans criterion. 
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4.2.3. Exploring Parameter Space 

Wc now harness the key advantage of our method, the 
abihty to perform controhed experiments, and ask how 
sensitive the results, obtained in Run A, are to variations 
in the initial conditions. Again, we here discuss the sim- 
ulations up to the formation of the first clumps, and turn 
to the further clump evolution later. 

(i) Spectral index 

To investigate the role of the DM substructure, we 
consider in Run E the case of a white-noise spectrum, 
P{k) oc k^. Admittedly, such a spectrum is physically 
ad hoc, in contrast to the k^^ case which is ultimately 
motivated by the theory of inflation. With the rms fluctu- 
ation on a given mass scale being (t(M) oc M~("+^^/^ for 
spectral index n, one finds 



a(M) cx 



const. for k ' 



(33) 
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Fig. 13. — - Run E: Morphology at z = 33.5. The convention in 
Fig. 7 is adopted for the rows and columns. Compared to Run A in 
Figure 8, there is relatively more DM substructure on the smallest 
resolvable scales, and the overall collapse proceeds in a more regu- 
lar way. The baryons do not yet fall into the shallow DM potential 
wells. 
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The k~^ case, therefore, has approximately equal power on 
all mass scales, whereas the smallest resolvable scale (given 
by the Nyquist frequency) dominates the white-noise real- 
ization. 

In Figures 13 and 14, we present the evolution of Run 
E up to the onset of gravitational instability. This time 
sequence is to be compared to the corresponding Figures 7 
- 9 for Run A. At redshift z = 33.5, briefly before virializa- 
tion, the baryons have not yet begun to fall into the shallow 
DM potential wells, in contrast to Run A. Towards the end 
of the free-fall phase, at 2; = 31.2, the gas has settled into 
a ring-like central configuration with a morphology which 
is somewhat more regular than in Run A. This distribu- 
tion derives from the rather homogeneous c;haracter of the 
overall collapse which deviates much less from spherical 
symmetry than the k~^ case. High-density clumps are 
formed again with initial masses close to ^ lO'^Afo, and 
the full difference to Run A becomes manifest only during 
the later evolutionary stages, as will be discussed below. 

The occurrence of the regular ring structure in Figure 
14 might be an artefact of the initial conditions in Run 
E which are highly symmetric. Such a ring-like configu- 
ration, however, is not a typical result. In general, the 
resulting morphology is a somewhat accidental feature of 
our results, and is not important for the main conclusions 
of this study. 
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Fig. 14. — Run E: Morphology at 2 = 31.2. The convention in 
Fig. 7 is adopted for the rows and columns. The boxsize is 30 pc. 
Shown is the situation briefly before the formation of the first clump 
(enclosed in the small box). Again, the dark matter is loosing its 
substructure in the process of virialization, whereas the gas has set- 
tled into a ring-like, central configuration which has a very knotty 
appeajrance. 



(ii) Random realization of DM fluctuations 

To address the problem of poor fc-space sampling for the 
longest wavelength modes, and the resulting morphologi- 
cal variations in the dark matter, we now compare Runs 
A and K, and the corresponding Figures 9 and 15. Run 
K has the same parameters as Run A, but a dark mat- 
ter component which is perturbed according to a different 
realization of the Gaussian random process. The appear- 
ance of the central gas configuration at z = 31.2 is quite 
different indeed, with a much more extended, filamentary 
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distribution of gas in Run K. Also, the first two clumps to 
form have masses close to ^ 500 A^q, compared to the one 
I5OOM0 clump in Run A. We will see below, however, 
that these two cases later on converge to a rather similar 
state, despite the differences during the early evolutionary 
stages. 
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(ill) Angular momentum 

In Runs C, A, and D with initial angular velocities of 
uj = 0.1, 0.2, and 0.4, respectively, wc investigate the in- 
fluence of angular momentum (or spin) on the evolution 
of the primordial gas. Increasing the amount of spin has 
two main effects. First, the moment of virialization and of 
the onset of gravitational instability is delayed, leading to 
the sequence of collapse redshifts: Zyir ~ 31.7, 31.2, and 
29.8 for Runs C, A, and D. Second, the gas is less centrally 
concentrated, resulting in a reduced rate for the merging 
of clumps, as will be discussed in the following section. 

The first effect can be understood in terms of a straight- 
forward modification of the analytical top-hat model, 
adding to it the presence of angular momentum. By con- 
sidering the energy balance at turnaround, an estimate 
for the turnaround radius can be obtained as follows (in 
dimensionless units where G = M = R = 1): 



Rta 3 ' 2 ' 



(34) 



where Ui and Hi are the initial angular velocity and Hub- 
ble parameter, respectively. With the virial radius given 
by Ryir — l/2Rta, the redshift of virialization is approxi- 
mately 



Fig. 15. — Run K: Case with different realization of P{k) oc k^^ 
noise. Shown is the inorpiiology at 2 = 31.2, briefly after the for- 
mation of the first two clumps with masses 360, and 520Mq {small 
box). The convention in Fig. 7 is adopted for the rows and columns. 
The boxsize is 60 pc. This case is to be compared to Run A in Figure 
9. Both cases have the same initial conditions with the exception of 
the random realization of the DM fluctuation field. It can be seen 
that the resulting morphology is very different here, with a much 
more extended central gas configuration. 
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Fig. 16. — Run H: Case of halo with total mass of 2 X lO^Mg. 
Shown is the morphology of the gas for two subsequent times. Top 
row: Gas distribution at 2 = 27.2, briefly after the formation of 
a clump with M ~ 400Mq. Bottom row: Gas distribution at 
z = 26.3. By now, the clump has grown in mass to M ~ ISOOMq. 
Left panels: Face-on view. Right panels: Edge-on view. The box 
has a linear size of 15 pc. In the case of this low-mass halo, where 
the requirement for efficient cooling, tcool < *// > is only marginally 
satisfied, one clump forms in the center of the DM potential, and 
reaches a final mass of ~ 2000Mq . 



1 + Zyir - {1 + Zyir,nr){^ - "^O-^X^) ■ (35) 

Choosing the redshift of virialization in the absence of ro- 
tation as Zyir^nr — 31.8, niccly reproduces the numerical 
results. The presence of angular momentum, therefore, 
delays the collapse by reducing the binding energy of the 
halo. 

(iv) Halo mass 

In Run H, we study the collapse of a less massive halo 
of total mass 2 x IO^Mq. This case only marginally sat- 
isfies the Rees-Ostriker criterion (see Section 2.3). The 
smaller halo mass translates into a lower virial tempera- 
ture, Tyir ^ 2000 K, which in turn leads to a reduced effi- 
ciency of II2 cooling. Consequently, the condition for the 
termination of the free-fall phase, tcool ~ iff, is reached 
already at lower densities (10^—10^ cm~^). The gas, there- 
fore, goes through a prolonged phase of quasi-hydrostatic 
contraction, and attains a roughly spherical configuration. 
In Figure 16, we show the central gas distribution at two 
successive redshifts. Only one clump forms in the center of 
the cloud with an initial mass of ~ 4OOM0. Subsequently, 
this clump grows in mass up to ^ 2000Mq. 

The main difference to Run A, besides the formation of 
only one clump instead of a few, is that here the Jeans 
instability proceeds less violent, since the opposing effect 
of pressure is non-negligible in this case. 

(v) Baryon fraction 

In a universe with a significant contribution to the crit- 
ical density in the form of vacuum energy, a larger frac- 
tion of the matter resides in the baryonic component. For 
Qm = 1 ^ f^A — 0.3, one has pB — 0.20pm, where pm 
is the total density of matter. In Run L, we consider a 
halo with 20% of the mass in gaseous form. From the 
outset, this case behaves very different from the simula- 
tions with a low baryon fraction. As can be seen in Figure 
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17, the gravitational instability is already triggered during 
the free-fall phase in the three dominating DAI eonden- 
sations. Since the thermal properties of the gas are not 
very different from those in Run A, and there is a four 
times larger amount of available gas, gravity more easily 
overwhelms thermal pressure. The peculiar character of 
Run L continues into the later evolutionary stages, where 
an altogether larger fraction of the gas is able to condense 
into high-density clumps. 
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Fig. 19. — Gas properties in simulations with different initial 
conditions I. Temperature vs. hydrogen number density. All runs 
are shown at very nearly the same instant, at z = 31.0. (a) Fiducial 
case (Run A): Halo of total mass 2x 10^ Mq, collapsing at Zyir = 30, 
initialized with P{k) oc k~^, and including HD cooling, (b) Vary- 
ing the power spectrum (Run E): Same as (a), but P{k) oc fc". (c) 
Varying the cooling (Run F) : Same as (a) , but no HD cooling, (d) 
Same as (a), but different realization of the random density field 
(Run K). 



Fig. 17. — Run L: Morphology of case with high at z = 31.7. 
The convention in Fig. 7 is adopted for the rows and columns. The 
boxsize is 90 pc. Shown is the situation briefly after the formation 
of the first clumps with masses of 750, 850, and 144OM0. The halo 
is still in its initial free-fall collapse, and the DM has not yet virial- 
ized. Already in this early dynamical stage, the gas in the deepest 
DM potential wells has undergone runaway collapse. This behav- 
ior is in marked contrast to the cases with a low baryon fraction 
(Qs = 0.05). 
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Fig. 18. — Run G: Case of top-hat with 

Zyir — 20. Shown is 

the morphology at 2 = 20.6, briefly after the first clump has formed 
with a mass of 92OM0 (in the region marked by the small box). 
The convention in Fig. 7 is adopted for the rows and columns. The 
boxsize is 90 pc. Compared to Run A (a top-hat of the same mass 
collapsing at Zyir — 30) in Figure 9, the resulting gas configuration 
is much more extended. 



10000 



^ 1000 



10000 



(a) 


(b) 

i" ■-;-i-Si.--:"-. 


^^^^^^^^^^^^ ' ' ' 1 ' 




(0 


' 1 ' ' ' 1 ' ' ' 1 ' ' ' 1 ' ' ' 1 ' ' ' 1 ' 







02468 -2 02468 
log n„ [cm"^] log n„ [cm"'] 



Fig. 20. — Gas properties in simulations with different initial con- 
ditions H. Temperature vs. hydrogen number density, (a) Varying 
the number of particles (Run B): Simulation with large number of 
SPH particles, shown si z = 31.0. For clarity and ease of compari- 
son, only every eighth particle is plotted, (b) Varying the angular 
momentum (Run C): Low-spin case, shown at z = 31.7. (c) Varying 
the halo mass (Run H); Less-massive halo, shown at z = 27.2. (d) 
Varying the collapse redshift (Run G): Halo with z^i^ ~ 20, shown 
at z = 20.6. The thermodynamic behavior, displayed in Figures 23 
and 24, is very robust under variation of the initial conditions, and 
in each case, the gas attains characteristic values of temperature and 
density close to T ~ 200 K and n ~ 10^ — 10** cm~^, respectively. 

(vi) Collapse redshift 

In Run G, we consider the collapse of a halo with total 
mass 2 x IO^Mq, collapsing at Zyir — 20. Run A and 
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G correspond to a ^ 3(t and ^ 2a peak in the Gaussian 
random field, respectively. We show the situation briefly 
after the formation of the first clump in Figure 18, which 
should be compared to the equivalent stage for Run A in 
Figure 9. The initial clump mass is again ^ IQ^Mq, and 
the main difference between the two cases lies in the much 
more extended morphology of the gas in Run G, with a 
linear size of ~ 40 pc compared to ~ 10 pc in Run A. 
The larger extension of the central gas configuration is 
simply due to the smaller binding energy of the halo in 
Run G, with almost the same amount of initial rotational 
energy as in Run A. Otherwise, the evolution of the two 
simulations is very similar. 

Summarizing the results from our exploratory survey 
up to the onset of gravitational instability, two impor- 
tant lessons can be learned. First, the morphology of the 
collapse varies significantly among the different cases, de- 
pending on the initial conditions. Second, the thermody- 
namic behavior of the priomordial gas is very similar for 
all the cases studied, despite the morphological diversity. 
To demonstrate this, we show in Figures 19 and 20 the 
location of the individual SPH particles in the tempera- 
ture vs. density plane for eight different runs. As can be 
seem, in each case, the gas does attain the characteristic 
values of the temperature and density, T ~ a few 100 K 
and n ~ 10^ — 10^ cm~^. In this preferred region of 
T — n space, the evolution of the system slows down, al- 
lowing to imprint the corresponding characteristic Jeans 
scale of Mj ~ lO^M© onto the gas. From examining Fig- 
ures 19 and 20, it is also evident that the gas particles 
begin their runaway collapse with these characteristic val- 
ues, to swiftly attain a temperature of T ^ 1000 K at the 
theshold density of n = 10* cm~^, at which point they are 
incorporated into a sink particle (i.e., a clump). 

We turn next to the discussion of how the high-density 
clumps subsequently evolve, and of the processes through 
which they grow in mass, the accretion of surrounding dif- 
fuse material and the merging of clumps. 

4.3. Later Evolution 

At the end of the free-fall phase, the primordial cloud 
has fragmented into clumps with initial masses of Mci ~ 
10^ — IO^Mq. These clumps are the basic elements in the 
building-up of a spectrum of masses through the processes 
of accretion and merging. To address the complex dynam- 
ics which is involved in the shaping of the mass function, 
it is essential to be able to follow the systems's evolution 
over a few dynamical timescales. The technique of cre- 
ating sink particles allows us to do so, which is another 
important advantage of our approach. 

In the following, we first discuss the physics of accretion 
and merging, and then the resulting distribution of clump 
masses. 

4.3.1. Accretion and Merging of Clumps 

The next question to ask is: What determines the frac- 
tion of gas which ends up in high- density clumps, and 
how does the difference in this quantity between the var- 
ious simulations arise? With Mb = IO^Mq being the 
total amount of gas, and t^yn ~ lO'' yr the initial dy- 
namical timescale of the DM halo, the overall rate of con- 
version between diffuse gas and clumps can be estimated 



to be Mconv — MB/tdyn — 10 Mq yr . The conver- 
sion rate, Mconv, sets an upper limit to the star formation 
rate (SFR). Feedback effects from the developing proto- 
stars are expected to limit the SFR to a somewhat smaller 
value which at present is not known with any certainty. 
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Fig. 21. — Growth of clumps in different runs. Solid line: Mass 
of most massive clump vs. time. Dashed line: Mass of second 
most massive clump vs. time. Mass is plotted in units of IO^Mq, 
and time in 2 X lO^yr, wrhich corresponds to the initial dynamical 
timescale. (a) Run A: Fiducial case, (b) Run C: Low-spin case, (c) 
Run E: P{k) oc fe". (d) Run K: Different realization of P(fc) oc k'^. 
A rapid rise in mass (a 'step') corresponds to a merger with another 
clump, whereas a phase of slow and steady growth is due to accre- 
tion of diffuse gas. Notice the similarity in panels (a) and (d), with 
two dominant clumps surviving without further merging. The sim- 
ulations in panels (b) and (c), corresponding to runs which result in 
a more centrally concentrated morphology, on the other hand, build 
up one very massive clump by successive merger events. 

The chimps are comprised of gas which at some point 
has become Jeans unstable, and which has been drawn 
from the reservoir of cooled gas at the characteristic val- 
ues of temperature and density, T ^ a few 100 K and 
n ~ 10^ — 10* cm^"^. In approximate pressure equilibrium 
with the cooled g second, hot phase has formed at 
T ~ 10** K and n ^ 10^ cm~^. Whatever material re- 
sides in this hot phase is not available for the formation 
of clumps. The total amount of cooled gas is typically 
Mcooi ~ 8 — 9 X 10'' M0. Once there is not enough gas 
left to overwhelm the opposing pressure, i.e., M < Mj, 
the Jeans instability ceases. If Mj is the average Jeans 
mass at the moment of virialization, where the average in- 
cludes all SPH particles which have been able to cool, the 
baryonic mass fraction in clumps is approximately given 
by _ 

Mcool - Mj 

In Table 4, we summarize the resulting conversion efii- 
ciency for three different simulations. The two simulations 

(Runs C and E) with a more concentrated gas morphology 
and correspondingly higher gas density are characterized 
by fci ~ 0-7, as compared to fci ~ 0.5 in the case of the 
less centrally concentrated gas configuration in Runs A 
and K. As can be seen, equation (36) nicely describes the 
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numerical results. All runs have a similar average temper- 
ature T ^ 300 K, but the average gas density is an order 
of magnitude larger in Runs C and E. For these latter 
two simulations, therefore, the corresponding Jeans mass 
is smaller, and less of the cooled gas is left behind in a pres- 
sure supported state. In general, the conversion efhcicncy 
fci increases with the central concentration and density of 
the gas. The larger central concentration is due to the low 
degree of angular momentum in Run C, and to the ab- 
sence of a significant deviation from spherical symmetry 
in Run E. The fraction fci constitutes an upper limit for 
the star formation efficiency (SFE), with the same degree 
of uncertainty as in the case of Mconv and the SFR. 

In Figure 21, we show how the two most massive clumps 
grow in mass for four different simulations, corresponding 
to the cases in Table 4. The clumps are formed with ini- 
tial masses close to M ~ lO^M©, and then gain in mass by 
the slow accretion of siirrounding gas, and by merging with 
other clumps. Both mechanisms can be clearly discerned 
in the figure, where the merging events correspond to the 
step-like, sudden increase in mass. It is evident that there 
is significant merging activity in the high-density simu- 
lations (Runs C and E), whereas in Runs A and K, the 
clumps at late times only grow by steady accretion. This 
difference can be understood by considering the timescale 
for the collision of clumps, tcoii, and the corresponding 
collision rate (e.g., Bonnell et al. 1998) 



the merging history of the clumps, is therefore expected to 
vary significantly. In the following, we investigate how the 
distribution of clump masses depends on the initial condi- 
tions. For each simulation, the gas and clump morphology 
is shown at a late evolutionary stage, where most of the 
available, cooled gas has already been incorporated into 
clumps. 
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Here, nci is the number density of clumps, v the velocity 
dispersion, race the accretion radius, and Adci the mass of 
a chimp. The second term in the brackets is the Safronov 
miml)er. describing the effect of gravitational focusing. 
Estimating the accretion radius as Tacc — '^hci — 0.1 
pc (see Section 3.2), and taking the clump mass to be 
Mci — 20, OOOMq, the result of evaluating equation (37) 
is summarized in Table 5. We accordingly expect of order 
5 merger events in the high-density simulations (Runs C 
and E), as opposed to only 1 event in the low-density cases 
(Runs A and K). As can be seen in Figure 21, this predic- 
tion is borne out in the numerical simulations. A clump 
can become very massive by undergoing multiple mergers, 
up to ~ 50, OOOMq in Run C, and ~ 60, OOOMq in Run 
E. This runaway growth of one central clump is analogous 
to the evolution of a supergiant cD galaxy in the center of 
rich clusters of galaxies. From Figure 21, it is evident that 
at late evolutionary stages, clumps grow in mass mainly 
by merging with other clumps. In comparison, accretion 
of surrounding gas is a rather inefficient process. 

4.3.2. Distribution of Clump Masses 

In this section, we discuss the clump mass spectrum 
which results from the complex dynamics of merging and 
accretion described above. Although the initial masses 
of the clumps are close to ~ IO^Mq, rather independent 
of the initial conditions, the subsequent evolution of the 
clumps proceeds differently from case to case. As we have 
seen, the efficiency of merging is very sensitive to the cen- 
tral density of the post-virialization cloud. Despite the 
existence of a characteristic mass scale for Population III 
star formation, the mass spectrum, being determined by 



Fig. 22.— Run C: Low-spin case at z = 28.9. Top row: The 
remaining gas in the diffuse phase. Bottom row: Distribution of 
clumps. The numbers next to the dots denote clump mass in units 
of Mq. Left panels: Face-on view. Right panels: Edge-on view. 
The length of the box is 30 pc. Dominated by a massive clump of 
~ 40, OOOMq, comprising ~40% of the initially present gas, a com- 
pact, disk-like feature has formed in the center of the DM potential. 
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Fig. 23. — Run A: Intermediate-spin case at z = 28.9. The 
manner of presentation is the same as in Figure 22. Compared to 
the low-spin case, the gas has settled into a less regular, more ex- 
tended configuration with two dominant clumps of mass close to 
20, OOOMq. During the subsequent evolution, the clumps survive 
without merging, and grow in mass only slightly by accretion of 
surrounding gas. 
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Fig. 24. — Run D: High-spin case at z = 25.5. The manner 
of presentation is the same as in Figure 22. Notice that here the 
situation is shown at a much later instant, reflecting the delayed 
conversion of the diffuse gas into clumps. The gas morphology is 
highly irregular and dispersed. In this case, no very massive clumps 
have formed, but instead a number of clumps with masses of a few 
10^ Mq. 

The total angular momentum of a halo proves to be a 
very important parameter in determining the mass spec- 
trum of the clumps. To demonstrate this, we consider the 
simulations with low, intermediate, and high initial spin 
in Figures 22, 23, and 24, respectively. As can be clearly 
seen, the mass of the dominating clump decreases with 
increasing spin. This clump mass - spin relation has a 
straightforward explanation in the lower central density 
of the higher spin simulations, and the corresponding re- 
duction in the merger rate. Hence, high-density clumps 
with masses M ^ 50, OOOM0 can form in the center of the 
lowest-spin halos. These massive clumps might conceiv- 
ably lead to the formation of seed black holes for future 
quasar activity. Eisenstein & Loeb (1995) have investi- 
gated a scenario along these lines, emphasizing the cosmo- 
logical importance of the low-spin systems. 

The emergence of such very massive (cD like) clumps 
is due to the efficient outward transport of angular mo- 
mentum via tidal torques. Although we do not include 
the effect of any external tidal fields, and have to insert 
the angular momentum explicitly at the beginning of the 
simulation, the subsequent dynamics does lead to internal 
tidal fields. In cases where there is such a " cD behavior" , 
the effect responsible for it therefore seems to be modelled 
in a physically consistent way. Interestingly, a high degree 
of central concentration is also found in the simulation of 
ABN2000 who have self-consistently included cosmologi- 
cal tidal fields. As we discuss below, the neglect of any 
proto-stellar feedback seems to be the more crucial factor 
in accounting for the emergence of very massive clumps. 

Earlier on, we have mentioned that the two simulations 
(Runs A and K) with a different realization of the 
power spectrum, and all other parameters being equal, 
lead to a rather different DM and gas morphology at the 
moment of virialization. During the later stages of the 
evolution, however, these simulations converge to a similar 
state. In both simulations, the morphology is dominated 



by two clumps with masses close to ~ 20,000AfQ. The 
overall features of the mass spectrum, therefore, do not 
seem to depend significantly on the random initialization 
of the DM fluctuations. The late-time morphology in the 
simulation with a high baryon fraction (Run L) is char- 
acterized by the largest fraction of gas in clumps, ~ 75 
%, and has a central cD clump of mass M ~ 200, OOOMq, 
which has gobbled up a remarkable one-half of the total 
gas mass. As we have pointed out before, the case of a 
high baryon fraction behaves distinctly different from the 
low l^B simulations in that gravitational instability is al- 
ready triggered during the initial free-fall phase. Now we 
see that Run L is also distinguished by an extremely top- 
heavy spectrum of clump masses. This case serves as an 
illustration for what happens when the relative balance 
of gas pressure and self-gravity is shifted in favor of the 
latter. 

To assess the possible role of HD cooling in the evolu- 
tion of the primordial gas, we compare Run A, which does 
include HD cooling, and Run F, which does not. Other- 
wise, the two cases have identical initial conditions. We 
find that the gas and clump morphologies, as well as the 
thermodynamic behavior in the T vs. n plane, are overall 
very similar. Slight differences, however, do exist in the 
way the gas fragments, and in the resulting clump masses. 
Cooling due to the HD molecule, and the corresponding 
chemistry of its formation and destruction, should there- 
fore be included for completeness, even though it does not 
appear to be of pre-eminent importance in the parameter 
regime considered here. In our simulations, we have as- 
sumed a deuterium abundance of no = 4 x 10~^nH, but we 
have carried out a test calculation with a ten times higher 
abundance. With such a high D abundance, the thermal 
evolution of the gas proceeds very differently. Cooling now 
is so efficient that the gas quickly settles down to the tem- 
perature of the CMB, and the corresponding Jeans mass 
is reduced below the resolution limit of that simulation, 
Mres ^ 200 A/q. If deuterium were indeed that abimdant, 
the characteristic mass scale of ~ IO'^Mq would then dis- 
appear, since it derives from the properties of H2 (see Sec- 
tion 2). Our adopted value of nD/nn = 4 x 10~^, however, 
seems to be a conservatively high choice, according to re- 
cent observations of high-redshift absorption systems (e.g., 
Buries k Tytler 1998). 

Finally, we discuss the cumulative distribution of clump 
masses. Wc consider the average of all simulations with 
DM fluctuations imprinted according to a spectrum, 
initial spin uj = 0.2, and a total halo mass of 2x IO^Mq. As 
we have discussed above, varying these parameters leads 
to a significant change in the distribution of clump masses, 
and we have therefore included in the averaging only the 
cases with these fiducial values. We evaluate the num- 
ber of clumps per unit logarithmic mass, and compare 
this to the Salpeter case of the present-day stellar IMF, 
diV/dlogm oc m~^'^^. It turns out that the resulting 
mass spectrum is remarkably flat, and that the halos in 
our simulations are very efficient in building up very mas- 
sive clumps, up to a few times lO^M©. 

This mass spectrum has to be taken cum grano salis, 
and important caveats apply. First, recall that clumps 
are born with masses close to ~ IO^Mq, and only sub- 
sequently reach higher masses through successive mergers 
and accretion. For these mechanisms to be effective, there 
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has to be sufficient time. Typical merging timescales are 
of order a few 10® yr, dming which time the first stars 
in the halo might have already evolved far enough to ex- 
plode as supernovae, thereby considerably disturbing the 
remaining gas in the halo. In general, the neglect of any 
feedback effects constitutes the major shortcoming in our 
treatment of these later evolutionary stages. Up to the 
point of forming the first clump, and possibly the first few, 
all the relevant physics is in principle known. It is only to- 
wards the later stages in the evolution that the physical 
basis of the simulations becomes increasingly uncertain. It 
is also worth remembering that only those processes are 
allowed for in our numerical treatment which tend to in- 
crease the mass of a clump. This introduces a bias towards 
higher mass in the resulting mass spectrum. 

Having stated all these provisions, it is nevertheless re- 
markable how efficient the Population III halos are in 
building up clumps of very high mass. The assessment 
of how relevant this finding is for an understanding of the 
stellar IMF, has to await a more realistic treatment of the 
complex physics of accretion from a dust- free envelope, the 
merging of subcondensations, and the various feedback ef- 
fects. 

5. SUMMARY AND CONCLUSIONS 

We have investigated the collapse and fragmentation of 
primordial, metal-free gas. The gas is embedded in dark 
matter halos of mass close to ^ IO^Mq which virialize 
at redshifts z ~ 20 — 30. Due to the low virial temper- 
atures in these halos of a few 1000 K, cooling can only 
proceed via H2. We find that for these systems there ex- 
ists a preferred region in parameter space. The primordial 
gas docs attain characteristic temperatures of a few 100 K, 
and densities of ~ 10^ — 10** cm~"^. These values have their 
physical explanation in the microphysics of H2 cooling, re- 
lated to the minimum temperature that can be reached 
via H2, and to the critical density where the cooling rate 
changes from being proportional to to an only linear 
dependence on density. This change in the cooling rate re- 
flects the transition from NLTE to LTE level populations 
of the hydrogen molecule. With these values of temper- 
ature and density, the corresponding characteristic Jeans 
mass is Mj ~ lO^M©. At some point during the simula- 
tion, the gas becomes gravitationally unstable, and forms 
high density (n > 10^ cm^^) clumps with initial masses 
close to the characteristic Jeans scale of ^ 10'^ M0. This 
is a very robust result, and is quite independent of the 
initial conditions. This suggests that Population III star 
formation might have favored massive stars, possibly even 
very massive ones with M ^ IOOMq. 

In contrast to this robust nature of the thermodynamics, 
the resulting morphology is a somewhat accidental feature 
of our simulations, varying substantially between different 
random realizations of the initial DM density field. 

Although the clumps form initially with similar masses 
close to the characteristic value of ^ 10'^ M©, their fur- 
ther evolution is very sensitive to the initial conditions. 
Clumps grow in mass through accretion of surrounding 
gas, and merging with other chimps. Both mechanisms 
sensitively depend on the central gas density. We now dis- 
cuss the dependence of the resulting clump masses on the 
most important parameters. 

Total halo mass. — With increasing mass of the DM halo, 



the clump masses also tend to increase. Initial clump 
masses arc close to Mci = 400il/,;) in the 10"'' M0 halo, 
and approximately twice as massive in the IO^Mq cases. 
In addition, a few clumps form almost simultaneously in 
the more massive halos, whereas only a single one forms 
in the center of the 'marginal' (10^ Mq) halo. 
Angular momentum. — The growth of a clump by accre- 
tion and merging is very sensitive to the amount of angular 
momentum initially imparted to the DM halo. Lower spin 
halos acquire denser central configurations which favor the 
frequent merging of clumps. Clump masses can then be- 
come very large (up to a few times 10* Mq). 
Collapse redshift. — The efficiency of accretion and merg- 
ing is also increased in runs with higher collapse redshifts, 
which is again a direct consequence of the enhanced over- 
all density. As a result, clump masses tend to be larger at 
higher redshift. 

Baryon fraction. — In the simulation with a larger fraction 
of baryons relative to the DM, runaway collapse is already 
triggered during the initial free-fall collapse, and resulting 
clump masses are systematically higher. 

Recently, ABN2000 have used the adaptive mesh refine- 
ment (AMR) method to study the formation of the first 
stars. These authors start their calculation with cosmolog- 
ical initial conditions, resulting in the most realistic treat- 
ment of the problem to date. Since they have presented 
only one case, however, it is not easy to ascertain the ro- 
bustness of their result, and this is the specific advantage 
of our exploratory approach. We agree on the existence 
of characteristic values for the density and temperature of 
the primordial gas, deriving from the microphysics of H2 
cooling. ABN2000 give an upper limit for the final mass of 
a Population III star of 2OOM0, whereas the distribution of 
clump masses in our work suggests that somewhat higher 
masses, perhaps up to IOOOMq, are possible. ABN2000 
have criticised our earlier work (Bromm et al. 1999) in 
that the rather smooth gaseous disk found in that study 
is unrealistic and an artefact of the assumed top-hat ini- 
tial conditions. Most of the cases presented here, however, 
do not yield such a regular disk morphology, and result in 
very irregular configurations that arc dominated by fila- 
ments and knots. Furthermore, our Run H, corresponding 
to a less massive halo {Mtot ~ lO^M©), is in good over- 
all agreement with the simulation of ABN2000 as to the 
morphology and the fact that only one clump forms in the 
center of the DM halo. 

Although one has a reasonably strong case in arguing 
that all the relevant physics is in hand to follow the col- 
lapse of the primordial gas up to the formation of the first 
high-density region, this fortunate circumstance gives way 
to growing uncertainty afterwards. Two of the most ten- 
talizing questions are the following. How effective is the 
energy output from the protostellar core in shutting off ac- 
cretion from a very massive, dust-free envelope? On this 
question hinges the ultimate mass scale of the first stars, 
and it is at present far from being answered. In Paper 
II, however, we will further address the question of the up- 
per mass limit of a Population III star. This upper mass 
limit is important for testing the idea that Population 
III stars might be the precursors of black holes of mass 
^ 10^ Mq, intermediate between stellar and supermassive 
ones (e.g., Madau & Rees 2001). What is the IMF of Pop- 
ulation III stars? In our simulations, we have seen how the 
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mass spectrum of clumps is built up through the complex 
interplay of merging and accretion, and that the resulting 
mass spectrum is very flat compared to the Salpeter IMF. 
It is difficult to tell, however, how this result relates to the 
stellar IMF, since again we are hampered by our neglect 
of any negative feedback effects, which are almost certain 
to play an important role in determining the IMF. Both 
questions highlight the importance of feedback effects for 
a deeper understanding of the star formation process. 
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APPENDIX 
A. N-BODY SOLVER 

In the following, we briefly discuss the basic ideas of solving the gravitational N-body problem, and refer the reader 
to Hernquist & Katz (1989; HK89 henceforth), Aarseth (1994), and Barnes (1998) for further details. In the context of 
TREESPH, the treatment of self-gravity is almost identical for the dark matter and the gas, since both components are 
represented by particles which constitute a Monte-Carlo sampling of the underlying fluids. Let us first turn to the dark 
matter, and describe the minor modifications for the gas at the end of this section. 

Since particle methods are Lagrangian by design, mass is conserved automatically, thus rendering the equation of 
continuity superfluous. The equation of motion for DM particle i simply reads: 

f = -(V*). . (Al) 

The gravitational potential is given by the standard solution to Poisson' s equation: 

$, = -G/^^dV . (A2) 

J \r-n\ 

If one now were to assume the case of true point masses, the density could be written as 

p{f^ = ^mjSir-fj) , (A3) 

j 

with 5{r} being the Dirac delta-function. Inserting equ. (A3) into equ. (A2), and applying the gradient operator, yields 
the familiar result for the gravitational acceleration of particle i: 

-(V$), = G^m,^f^ . (A4) 

This straightforward procedure has two crucial shortcomings which we now discuss in turn, together with the adopted 
remedies. 

First, there is the problem of discreteness on very small scales. Close encounters of particles lead to collisional two-body 
relaxation, which is clearly undesirable in modelling a smooth fluid. Therefore, gravitational forces have to be softened on 
scales close to the interparticle distance. In TREESPH, this is accomplished by replacing the singular Dirac delta-function 
with the spherical spline kernel W, used in the SPH formalism (see Appendix B): 

p{f) = Y^mjW{f-fj;e) . (A5) 
j 

The result of inserting this Ansatz into equ. (A2) is given by HK89. Due to its compact nature (W(r) = for r > 2e), 
the use of this kernel reproduces the Kepler-potential exactly for r > 2e. Each particle has its own softening length, Cj, 
which is adjusted such that there is a constant number of neighbors, Ns, in a softening volume: Ns rij (2ej) , where rij 
is the local particle density. In our simulations, we have set Ng = 16. 

The second problem concerns the computational expense of determining the gravitational forces, which scales as 0{N'^) 
for the direct summation in equ. (A4). An ingenious way to overcome this prohibitive cost has been developed by Barnes & 
Hut (1986). Their method organizes the particles into a recursive tree structure, such that individual particles correspond 
to the leaves of the tree, neighboring particles in space to leaves on the same branch, and the system as a whole to the 
root of the tree. Each node of this tree represents a cubical cell in real space. The leaves correspond to the smallest cells, 
containing only one particle, higher-level branching points to larger cells, containing groups of particles, and the root to 
the overall system box, including all the particles. This procedure is completely adaptive, and can accomodate arbitrary 
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geometries. To finally calculate the gravitational force on a given particle, the tree is traversed from the root down. For 
each node (cell) of size s and distance d to the particle, an accuracy criterion is evaluated: 

, (A6) 

where the parameter 9 determines the desired precision, and is chosen to be = 0.8 in our simulations. If this criterion 
is fulfilled, all the particles in the given cell are treated as a single group whose gravitational influence is approximately 
described by performing a multipolc expansion of the respective potential, in our case up to quadrupolc order. If, 
on the other hand, the accuracy criterion is not met, the tree is walked down to the next level of refinement. This 
procedure is repeated recursively, until either the criterion in equ. (A6) is fulfilled, or the tree descent has reached the 
level of individual particles (the leaves) . In consequence, forces are calculated accurately but expensively only for nearby 
particles, and approximately but fast for more remote ones. Both the construction of the tree, and the subsequent descent 
of it, have a computational cost of only 0{NlogN), making simulations with larger numbers of particles possible. 

The equation of motion is solved explicitly using a standard, time-centered leapfrog integrator. The position and 
velocity of particle i are updated according to: 

^+1/2 ^ -^-1/2 ^ ^^^^ ^ Q ^^^3) 

tJ^+i =v^ + a"+^/^At, + O (Atf) (A7) 

Here, superscripts denote timesteps. Velocities are stored at times which are offset by one-half a timestep from the 
positions and accelerations. This (leapfrog) characteristic guarantees the second-order accuracy of the algorithm. Each 
particle is allowed to have its individual timestep, Afj, which is chosen to fulfill the criterion 

Ati < etoi— . (A8) 
aiVi 

Ei is the total energy of particle z, and etoi the tolerance parameter which we have set to be ctoi = 0.1. 

The gaseous component is treated in exactly the same way, with the exception that in this case the gravitational 
softening length is always equal to the SPH smoothing length (see below). In Figure 1, we demonstrate that the code 
does nicely reproduce the analytical top-hat solution, with total energy and angular momentum being conserved to better 
than 5%. 



B. THE SPH METHOD 

We briefly describe the basic principles of the smoothed particle hydrodynamics (SPH) method, and give an intuitive 
motivation for the resulting equations. Further details are given in, e.g., Benz (1990), Monaghan (1992), Miiller (1998), 
and HK89. 

The fluid is sampled by discrete particles, representing fluid elements. To model a continuous medium, the mass of a 
particle is smoothed according to 

p{f) = Y,mjW{r-rf,h) (Bl) 
j 

The smoothing kernel, W , is strongly peaked at r = r^, and normalized to give J W^d^r = 1, where integration is over all 
space. TREESPH implements the spherically symmetric spline kernel 

1 f + o<^<i 

^^'^''^ = ^^\ 1<^<2 (B2) 

I 0, i>2 

The smoothing length, h, describes the spatial extent of a given SPH particle. TREESPH assigns variable smoothing 
lengths to each particle, thereby introducing an adaptive spatial resolution, such that there is a (roughly) constant number 

of particles, Nndgh, within the smoothing volume. We have adopted Nndgh = 32 in our work. Now, one can easily find 
smoothed estimates for any variable A{r). Starting with the interpolation formula 

A{r) = j A{r')W{r-r';h)d?r' , (B3) 

and assigning an effective volume to each particle j via rrij ~ p{fj)A^rj, one can approximate the integral as 

A(r)^5]A,^I^(f-r,;/i) , (B4) 
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where Aj = A{rj), and pj = p{rj). The crucial advantage of the SPH method Ues in the fact that it does not need a grid 
to evaluate spatial derivatives. To illustrate this point, let us apply the gradient operator to equ. (B4): 

VrA{f) ~ ^ Aj-^VrW(f - fjih) . (B5) 
J 

Consequently, taking spatial derivatives in SPH only involves the analytical differentiation of the kernel function which is 
specifically chosen to admit this. 

Next, we present the SPH equations as used in our numerical work. We give heuristic arguments for their appearance, 

and refer to the cited literature for the formal derivations. The equation of continuity is again automatically fulfilled due 
to the Lagrangian nature of the SPH method. The equation of motion for particle i reads 

dvi VP,- 

d^ = -(^*)^-7f + "^-'' • ^^'^ 

To find a smoothed estimate for the pressure gradient, one could naively write, using equ. (B5): 

VP- P 

Pi y P^P3 

This form, however, does not conserve linear and angular momentum. To satisfy Newton' s third law, one wants an 
expression which is symmetric with respect to any given pair of particles, i and j. Prom amongst the options given in 
TREESPH, we have chosen to work with 

VP / p p 1 

^ = E 9 ["^iWinfM + y^w{nf,hi)] . (B7) 

Pi ^. PiPj ^ 

The symbol Vj denotes derivation with respect to ri, and rij = \ri — fj\. The pressure is given by the ideal gas law 

Pi = = {j - l)uipi . (B8) 

IJ.mil 

jj, is the (dimensionles) mean molecular weight, 7 the ratio of specific heats, and the specific internal energy (in erg g~^). 

For an almost neutral gas consisting of helium and atomic hydrogen, one has p ~ 1.2, and 7 = 5/3. The corresponding 
symmetric estimate for the viscous acceleration dsh = — VQeff/p can be written as: 

= V Ilijmj I [ViWinj ; hi) + ViWinj ; hj)] . (B9) 
Pi . I 

Finding Hjj, a symmetric estimate for Q^nj ^ relies on the following conceptual steps. First, let us specify the pseudo- 
viscous pressure, QeS = —p'^es'^ ■ arising from the presence of an artificial viscosity. The latter has two contributions: 
A bulk viscosity 

Ves = alcs , (BIO) 

and a von Neumann-Richtmyer viscosity 

u^B=-pl^W-v . (BU) 

Here, I is a characteristic length, Cg the sound speed, and a, f3 are parameters of order unity. The resulting pseudo- viscous 
pressure can then be written as 

QeB = -apCslV-v + (3pl''{V -vf . (B12) 
A symmetric estimate for IV ■ v is given by 

= } FpOiXf- ^ u ^ (B13) 

[ otherwise 

where Vij = Vi — Vj, and hij = {hi + hj)/2. Assembling these ingredients, one finally has 

-api,ci, + Pp^, 

Pij 

where pij = {pi + pj)/2, and c^- = {cs^i + Cs,j)/2. In our work, we use a = 1 and (3 = 2. We have experimented with 
alternative prescriptions for the artificial viscosity (as given by HK89), and find that our results do not depend on this 
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choice. Certain morphological features such as pronounced rings, however, could be due to an unphysical viscosity that 
can arise in simulations with a relatively small number of particles. The equation of motion is now fully specified, and is 
solved with the leapfrog integrator. The stability of the integration is enforced by the classic Courant condition: 

^U< • (B15) 

CsA + rii\\/ ■ Vi\ 

HK89 give a version of the criterion which also takes into account the artificial viscosity. The appropriate timestep for a 
given SPH particle is the minimum of the Courant timestep and the one given by the energy criterion, equ. (A8). 
Along similar lines, HK89 write the smoothed form of the thermal energy equation, equ. (5), as 



■l[ViWirif,hi)+\/iW{ri,;h,)] 



r- A 



In the presence of radiative cooling and heating, this equation cannot be solved explicitly, since the corresponding radiative 
timescalcs are typically much shorter than the dynamical time which sets the timestep for the equation of motion. The 
thermal energy equation is therefore integrated implicitly, using the standard second-order trapezoidal rule (see HK89 for 
further details). 
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Table 1 

Reaction Rates for Hydrogen Species 







Reaction 




Rate Coefficient 
(cm^s"!) 


Reference 


(1) 


H 4 


- e- ^ H+ 4- 


■ 2e- 


5.85 X 10-"Ti/2cxp(-157, 809.1/T)(1 +T5/^)-i 


1 


(2) 


H+ 


+ e- H + 


hu 


8.40 X 10""T-1''2j.~o.2j.j^ _j_y0.7-)-i 


1 


(3) 


H ^ 


- e- + hu 


See expression in reference 


2 


(4) 


H 4 


- H- ^ H2 + 


- e 


1.30 X 10^3 


1 


(5) 


H- 


+ H+ ^ 2H 




7.00 X lO-^T-i/2 


1 


(6) 


H2 


+ e- ^ H + 


R- 


2.70 X 10"*^T-3/2cxp(-43,000/T) 


1 


(7) 


H2 


+ H ^ 3H 




See expression in reference 


1 


(8) 


H2 


+ H+ ^ H+ 


+ H 


2.40 X 10--'cxp(-21,200/T) 


1 


(9) 


H2 


+ e- ^ 2H H 


h e- 


4.38 X 10""'cxp(-102,000/r)rO-38 


1 


(10) 


H- 


+ e- ^ H + 


■ 2e- 


4.00 X 10-i2Texp(-8750/T) 


1 


(11) 


H- 


+ H ^ 2H 4 


- e~ 


5.30 X 10-20Texp(-8750/T) 


1 


(12) 


H- 


+ H+ ^ H+ 


+ e- 


See expression in reference 


1 



References. — (l)Haiman, Thoul, & Loeb 1996; (2) Abel, Anninos, Zhang, & Norman 1997. 
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Table 2 

Reaction Rates for Deuterium Species 





Reaction 




Rate Coefficient 

(cm^s^-*-) 


Reference 


(1) 


D+ + e- ^ D + hi/ 


8.40 X 10- 


-llj,~l/2-p-0.2(^_^y0.7)-l 


1 


(2) 


D + H+ ^ D+ + H 


3.70 X IQ- 


-l0T0-28exp{-43/r) 


3 


(3) 


D+ + H ^ D + H+ 


3.70 X IQ- 


-102^0.28 


3 


(4) 


D+ + H2 ^ H+ + HD 


2.10 X IQ- 


-9 


3 


(5) 


HD + H+ ^ H2 + D+ 


1.00 X IQ- 


-''exp(-464/T) 


3 



References.— (l)Haiman, Thoul, & Loeb 1996; (3) Galli & Palla 1998. 



Table 3 

Parameters for the Different Runs 





M/Mq 




n 


U! 




HD 




Run A 


2 X 10^ 


30 


-3 


0.2 


0.05 


Yes 


16384 


Run B 


2 X 10^ 


30 


-3 


0.2 


0.05 


Yes 


131072 


Run C 


2 X 10« 


30 


-3 


0.1 


0.05 


Yes 


16384 


Run D 


2 X 10^ 


30 


-3 


0.4 


0.05 


Yes 


16384 


Run E 


2 X 10« 


30 





0.2 


0.05 


Yes 


16384 


Run F 


2 X 10« 


30 


-3 


0.2 


0.05 


No 


16384 


Run G 


2 X 10^ 


20 


-3 


0.2 


0.05 


Yes 


16384 


Run H 


2 X 10^ 


30 


-3 


0.2 


0.05 


Yes 


16384 


Run K" 


2 X 10^ 


30 


-3 


0.2 


0.05 


Yes 


16384 


Run L 


2 X 10^ 


30 


-3 


0.2 


0.20 


Yes 


16384 



"Run K has the same parameters as Run A, but with a diflferent reaUzation of the random density 
field. 



Note. — n is the spectral index, u) the dimensionless angular velocity, Qb the baryon fraction, and 
HD refers to the absence or presence of HD cooling. 



Table 4 

Efficiency of Forming Clumps 





^dyn 




<Mj> 


fci 




[yr] 


[Mq] 


[Mq] 




Run A 


6.0 X 10*5 


8.0 X 10-* 


3.0 X Iff* 


0.5 


Run C 


5.4 X lO^ 


9.0 X 10-* 


2.0 X Iff* 


0.7 


Run E 


3.0 X 10" 


8.5 X 10* 


1.5 X Iff* 


0.7 



Note. — Parameters for the simulations in Figure 21. tdyn is the dynamical time of the DM halo 
at the moment of virialization, M^^gi the total amount of gas which has been able to cool, < Mj > 
the average Jeans mass at virialization, and fci the baryonic mass fraction in clumps. Run K has 
virtually the same parameters as Run A. 
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Table 5 
Expected Merging of Clumps 





nci 




tcoll 


^merger 






[km s^^] 


[y] 




Run A 


5 X 10-3 


12.3 


1.5 X 10'' 


~ 1 


Run C 


3 X 10-2 


10.3 


3.0 X 10® 


~ 5 


Run E 


3 X 10-2 


14.3 


3.0 X 10® 


~ 5 



Note. — Explaining the merging histories in Figure 21. nci is the number density of clumps, v 
the velocity dispersion, tcoii the collision timescale, and Nmerger the expected number of mergers in 
At ~ 2 X lO'^ yr. Run K has virtually the same parameters as Run A. 



